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ABSTRACT 

We use the results of magnetohydrodynamic (MHD) simulations to emulate spectroscopic observations and 
use maps of centroids to study their statistics. In order to assess under which circumstances the scaling prop- 
erties of the velocity field can be retrieved from velocity centroids, we compare two point statistics (structure 
functions and power-spectra) of veloc ity centroids with those of t he underlying velocity field and analytic pre- 
dictions presented in a previous paper (iLazarian & Esquive 1120031) . We tested a criterion for recovering velocity 
spectral index from velocity centroids derived in our previous work, and propose an approximation of the early 
criterion using only the variances of "unnormalized" velocity centroids and column density maps. It was found 
that both criteria are necessary, however not sufficient to determine if the centroids recover velocity statistics. 
Both criteria are well fulfilled for subsonic turbulence. We find that for supersonic turbulence with sonic Mach 
numbers Ais ^ 2.5 centroids fail to trace the spectral index of velocity. Asymptotically, however, we claim that 
recovery of velocity statistics is always possible provided that the density spectrum is steep and the observed 
inertial range is sufficiently extended. In addition, we show that velocity centroids are useful for anisotropy 
studies and determining the direction of magnetic field, even if the turbulence is highly supersonic, but only 
if it is sub-Alfvenic. This provides a tool for mapping the magnetic field direction, and testing whether the 
turbulence is sub-Alfvenic or super- Alfvenic. 

Subject headings: ISM: general — ISM: structure — MHD — radio lines: ISM — turbulence 



1. INTRODUCTION 

It is well established that the interstellar medium (ISM) 
is turbulent. From the theoretical point of view this arises 
from the very large Reynolds numbers present in the ISM 
(the Reynolds number is defined as the inverse ratio of the 
dynamical timescale -or eddy turnover time- to the viscous 
damping timescale). From an observational standpoint there 
is also plenty of evidence that supports a turbulent ISM, 
where the turbulence expands ov er scales that range from 
Au to kpc (Larson .1992: Arms trong. Rickett & Sr)an^ ien 
19951 iDeshpande Dwarakan ath. & Goss 20C3 
Stanimrrovic & Lazarianri200lTK This turbulence is very 



important for many physical processes, including star for- 
mation, cosmic ray propagati on, heat transport, a nd heating 
of the ISM (for a review see Vaz auez-Semadeni et aL.2000t 
ICho & Laza rian 2005; and references therein). 

How to compare interstellar turbulence with the results of 
numerical simulations and theoretical expectations is an im- 
portant question that must be addressed. After all, theoretical 
constructions involve necessary simplifications, while numer- 
ical simulations of turbulence involve Reynolds, and magnetic 
Reynolds numbers that are very different from those in the 
ISM. Are numerical simulations of ISM of value? To what 
extent do they reproduce interstellar turbulence? These sort 
of questions one attempts to answer with observations. 

Substantial advances in understanding of scaling of com- 
pressible MHD turbulence' (see reviews bv ICho & LazarianI 
|2()05; Lazarian & Cho 2005, and references therein) allow to 
provide a direct comparison of the theoretical expectations 
with observations. How reliable are the turbulence spectra 
obtained via observations? 

Studies of statistics of turbulence have been fruitful us- 
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' This scaling was used to solve important astrophysical probl ems, for 
instance, finding the rates of scat tering of cosmic ray s ( Yan & Lazarianl20 
and acceleration of cosmic dust lYan. Lazarian. & D raine 200J 



ing interstellar scintillations ("Naravan & Goodman' '19851 
ISpangler & Gwinn 1990) . However this technique is re- 
stricted to the study of ionized media, and very importantly 
to density fluctuations alone (see ICordeslll999h . Nowadays, 
radio spectroscopic observations of neutral media provide 
us with an enormous amount of data containing information 
about interstellar turbulence, including a more direct physi- 
cal quantity to study turbulence: velocity. But the emissivity 
of a spectral line depends on both the velocity and density 
fields simultaneously, and the separation of their individual 
contribution is not trivial. Much effort has been put into this 
difficult task and several statistical measures have been pro- 
posed to extract in formation of ve locity from spectroscopic 
data (see review bv lLazari'anlll999h . Among the simplest we 
can mention the use of line- widths |l arson 1981, 1992;Scal3 
11984, 1987"). Velocity centroids have been around as mea- 
sure of the velocity field for a long time now (^von HoerneJ 
1951; Munch 1958). And they hav e been widely used to study 
turbulence in molecular clouds ('Kleiner & Dick manI Il985t 
Dickman & Kleiner 1985; Miesch & Ballv 1994). 

Power-spectra, correlation, and structure functions have 
been traditionally, and still are, the most widely used tools 
to characterize the statistics of emissivity maps. These sta- 
tistical tools have been used to study the scaling properties 
of turbulence, e.g. to determine its the spectral index. Re- 
cently, more elaborated techniques have been proposed to 
analyze observational data, such as "A-variance" wavelet 
transfo rm (Stutzki et al. 1998; Mac Low & Ossenkopf 200(i 
lOssenkopf & Mac Low . 2002.) . Such techniques can be used 
to obtain velocity information from spectral data with some 
advantages, like being less sensitive to the effects of edges or 
noise. Regardless of the method used, it is usually assumed 
that the map traces the velocity fluctuations, which as we 
show below this is not always true. To separate velocity from 
density contribu tion "Modified Velocity Ce ntroids" (MVCs) 
were derived in (ILazarian & Esquivel2003[ hereafter LE03). 
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There has been an effort in parallel to develop new statis- 
tics that trace velocity fluctuations. Here we can mention 
the "Spectral Correlation Function" -SCF- fRosolowskv et al' 
[1999; Padoan, Rosolowskv & Goodman 2001), "Velocity 
Channel A nalysis" -VCA- (Lazarian & Pogosvan 2000; 
iLazarian et al . 200k lEsqu ivel et al. 2003), MVCs (LE03), 
and "Velocity Coordinate Spectrum" -VCS- (Lazarian 2004). 
Both VCA and VCS are good for studies of supersonic tur- 
bulence^. Although SCF was introduced as an empirical tool, 
its properties, in what the statistics of turbulence is concerned, 
can be derived using a general theory in L azarian & Pogosvan 
(|2004). 

Synergy of different techniques is very ad- 
vantageous for studies of inter s tellar turbulence. 
iMiville -Deschenes. Levrier & F algarond ( l2003bft attempted 
to test the results obtained with VCA using velocity centroids. 
However, as we show in the paper, without a reliable criterion 
of whether velocity centroids reflect the velocity statistics 
such studies deliver a rather limited insight. 

Numerical simulations provide us with an ideal testing 
ground for the statistical tools available for application to ob- 
servational data. However we must note that the situation is 
rather complex. On one hand, real observations depend criti- 
cally on the physical properties of the object under study, such 
as variations in the excitation state of the tracer and the radi- 
ation transfer within it (see Lazarian & Pogosvan 2004). In 
addition observational limitations, like finite signal-to-noise 
ratio and map size, griding effects, beam pattern, beam error, 
etc. are also present. On the other hand, numerical simula- 
tions have their own limitations, such as finite box size and 
resolution, numerical viscosity, and the physics available to a 
particular code. This paper is mostly concerned about the pro- 
jection effects and the impact of density fluctuations to cen- 
troid maps, which are shared by observations and simulations. 

In LE03 we studied the maps of velocity centroids as trac- 
ers of the turbulent velocity statistics. We derived analytical 
relations between the two point statistics of velocity centroids 
and those of the underlying velocity field. We also identi- 
fied an important term in the structure function of centroids 
which includes information of density, and that can be ex- 
tracted from observables. Subtraction of that term can isolate 
better the velocity contribution, and this yielded to a new mea- 
sure that we termed "modified velocity centroids" (MVCs). In 
LE03 we proposed a criterion for determining whether veloc- 
ity centroids reflect the scaling properties of underlying turbu- 
lent velocity (e.g. structure functions or spectra of velocity). 
A major goal of this paper is to test the predictions in LE03 
using synthetic maps obtained via MHD simulations and to 
determine when velocity centroids indeed reflect the velocity 
statistics. 

Earlier on, in ILazarian. Pogosvan & Esquivell (|2QQ3) we 
showed how velocity centroids can be used to reveal the 
anisotropy of MHD turbulence and how this anisotropy can 
be used for studies of plane-of-s ky magnetic field. This 
technique was further discussed bv Vestuto. Ostri ker & Stona 
(^003). In this paper we show how Mach number and Alfven 
Mach number affect the anisotropy of velocity centroid statis- 
tics. 

The results of LE03 obtained in terms of structure functions 
are trivially recasted in terms of spectra and correlation func- 
tions. Therefore we use structure, correlation functions and 

^ Using species heavier than hydrogen one can study subsonic turbulence 
as well 



spectra interchangeably through our paper, depending what 
measure is more convenient. While being interchangeable, 
for practical statistical data handling different measures have 
their own advantages and disadvantages. We discuss those 
on the example of power-law scalar field and thus benchmark 
our further velocity centroid study. We also deal with a po- 
tentially pernicious issue of non-uniformity of notations and 
normalizations that plague the relevant literature by having 
details of our derivations in the appendixes that constitute an 
important part of the paper 

In this work we perform a detailed numerical study of the 
ability of velocity centroids to extract turbulent velocity statis- 
tics. We study the issues of velocity-density correlations and 
outline the relation of velocity centroids to other techniques. 
In §2 we review the basic problem of the density and velocity 
contributions to spectroscopic observations. We summarize 
LE03 in §3, we include in this work appendices with math- 
ematical derivations omitted in our earlier short communica- 
tion. In §4 we test the analytical predictions, and the spectral 
indices from our numerical data. In §5 we show how cen- 
troids can be used for turbulence anisotropy studies and de- 
termination of the plane-of-sky direction of magnetic field. A 
discussion of the results can be found in §6, and a summary 
in §7. 

2. TURBULENCE STATISTICS AND SPECTRAL LINE 

DATA 

Due to the stochastic nature of turbulence it is best de- 
scribed by statistical measures. Among these we have two 
point statistics such as structure functions, correlation func- 
tions, and power spectra (see for instance Monin & Yagloi^ 
119751) . Their definition and a more comprehensive discussion 
can be found in Appendix A. Structure and correlation func- 
tions depend in general on a "lag" r, the separation between 
two points xi and X2, such that r = X2-X1. Power spectrum is 
defined as the Fourier Transform of the correlation function, 
and is function of the wave-number vector k. With ampli- 
tude ^ = |k| ^ 27r/r, where r = |r|. Additional simplification is 
achieved if the turbulent field is isotropic, in which case struc- 
ture and correlation functions depend only on the magnitude 
of the separation r (and not on the direction), similarly power 
spectrum is only function of k. This is not strictly true for 
magnetized media, as the presence of a magnetic field intro- 
duces a preferential direction for motion. In fact, MHD turbu- 
lence becomes axisymmetric in a system of reference defined 
by the direction of the local magnetic field (see reviews by 
Cho & Lazarian 2005 and Ch o. Lazarian & Vishn iac 2003), 
thus breaking the isotropy. However, since the local direction 
changes from one place to another, the anisotropy is rather 
modest and it is still possible to characterize the turbulence 
with isotropic statistics (see Esquivel et al. 2003). 

2.1. Three-dimensional power-law statistics 

In the simplest realization of turbulence we have injection 
of energy at the largest scales. The energy cascades down 
without losses to the small scales, at which viscous forces be- 
come important and turbulence is dissipated. At intermediate 
scales, between the injection and the dissipation scales, the 
turbulent cascade is self -similar. This range constitutes the 
so called inertial- range. There the physical variables are pro- 
portional to simple powers of eddy sizes, and the two point 
statistics can be desc ribed by power-laws. For power-law 
statistic s iLazarian fe'Pogosvan (2000) discussed two regimes, 
a short-wave-dominated regime, corresponding to a shallow 
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spectrum, and a long-wave-dominated specivum, with a steep 
spectrum. 

While dealing with numerical data one encounters a few 
non-trivial effects that we find advantageous to discuss below. 
The insight into the limitations of numerical procedures that 
involve conversion from mathematically equivalent statistics 
helps us for the rest of the paper 

2.1.1. Steep (long-wave-dominated) spectrum 

Consider an isotropic power-law one-dimensional power- 
spectrum of the form 



Lona — wave-doninated 



Pw(k) = C P" 



(1) 



A steep spectrum corresponds to spectral indices < -1. 
The structure function D{r) can be written in terms of the 
spectrum as 



D(r) = 4 / [1 -cos(;t r)] Pw(k) dk. 



(2) 



For a power-law steep power spectrum (substituting eq.fll 
into eq.|2l), the structure function also follows a power-law: 



D(r) = A r« = A r'^"''" < ^ < 2, 



(3) 



where A = C 27r/[r(l +7)sin(7r7/2)], and T(x) is the Euler 
Gamma function. The relation between the spectral index of 
the structure function and power spectrum can be generalized 
for isotropic fields to P^o oc k'"^", with 



1ND = -N-C 



(4) 



Where N is the number of dimensions (see Appendix A). For 
instance, the velocity (v) in Kolmogorov turbulence scales 
as V cx r'/-^, which corresponds to a spectral index for the 
structure function of ^ = 2/3, to a three-dimensional power- 
spectrum (Pio) index jjd = -11/3, a two-dimensional power- 
spectrum {P2d) index 72D = -8/3, and a one-dimensional 
power-spectrum Pio index jid = -5/3. Note that Kolmogorov 
spectrum falls into the steep spectra category. 

Structure functions given by equation (|3j are well defined 
only for ^ > 0, which allows to satisfy D(0) = 0, and f < 2, so 
that the representation in terms of Fourier integrals is possible 
(see Monin & Yaglom 1975).^ 

To illustrate the relation between the two point statistics and 
power-law spectrum, as well as the difficulties associated with 
handling numerical data. We produced a three-dimensional 
Gaussian cube with a prescrib ed (3D) spectral index of 73/5 = 
-11/3, as described in Esauivel et alJ (120031) . This type of 
data-cubes are somewhat similar to the fractiona l Br ownian 
motion (fBms) fields us ed by Brunt & Hever ( 20 02al) . or the 
de-phased fields used in Brunt et alj (i2003i) . However, as in 
real observations, they do not have perfect power-law spec- 
trum for a particula r realization, but only in a statistical sense 
(see .Esquivel et all 12003..) . In Figure [0 we show the calcu- 
lated 3D power spectrum, structure and correlation functions 
of our steep Gaussian cube, and compare them with the pre- 
scribed scaling properties. The power spectrum is computed 
using Fast Fourier Transform (FFT) in 3D and then averaged 
in wave-number k. Ideally one can compute directly in real 
space the 3D structure or correlation functions, however for a 

' Correlation functions in tliis regime are maximal, and finite at r = 0. It 
is important to notice also that for a power-law steep spectrum th e str ucture 
function is a power-law, but the correlation function is not (see eq. IA4I ). The 
correlation function is a constant minus a growing (positive index) power-law, 
therefore in a log-log scale is flat at small scales, and drops at large scales. 
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Fig. 1 . — Three-dimensional two point statistics of a long-wave dominated 
Gaussian field (steep spectrum). In panel (a) the three-dimensional power 
spectrum, the solid line is the prescribed spectrum, with a spectral index of 
730 =—11/3. the stars correspond to the calculated. In panel {b) the expected 
structure function (solid line), the computed structure function (stars); the 
expected correlation function of fluctuations (dotted line) , and the calculated 
correlation function of fluctuations(£(/amon&). 



3D field of the dimensions used here (216^) is already quite 
expensive computationally. It would require looping over all 
the points in the data-cube to do the average and repeat for all 
the possible values for the lag (in three dimensions as well). 
Fortunately, since the data-cubes were produced using FFT 
we can safely compute the correlation function with spectral 
methods. The correlation function can be expressed as a con- 
volution integral B{r) oc J dr f{r)f{r+r ), which can be cal- 
culated as a simple product of the Fourier tranformed fields. 
That is, B{r) oc T{f{k) f*(k)}, where /(k) = T{f(r)} is the 
Fourier transform of f{r), and f*ik ) its complex conjugate. 
Then, with the use of equation iA4i we can obtain the struc- 
ture function. The resulting correlation and structure func- 
tions in 3D are then averaged in r. An important thing to 
notice is that the Gaussian cubes have wrap-around period- 
icity, and the largest variation available correspond to scales 
of L/2, where L is the size of the computational box. In fact, 
we only plot the structure and correlation functions up to such 
separations. We see a fair agreement with prescribed and the 
measured scaling properties. The power spectrum in Figure 
nifl) shows departures from strict power-law which are more 
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evident for small wave-numbers. This is natural for this type 
of data-cubes, where random deviations from strict power-law 
are expected at all scales. But at large scales (small k) we have 
fewer points for the statistics and the departures do not aver- 
age to zero, while at small scales they almost do. 

2.1.2. Shallow ( short-wave-dominated) spectrum 

When the energy spectrum is shallow (i.e. 710 > -1), the 
fluctuations of the field are dominated by small-scales, there- 
fore termed short-wave-dominated regime. Density at high 
Mach number is an example of such shallow spectrum (see 
[Beresnvak, Lazarian & Cho 2005). In this case neither the 
structure nor correlation functions can be strictly represented 
by power-laws. In fact, in order for the Fourier Transforms to 
converge in this case we need to introduce a cutoff for small 
wave-numbers, such that the power spectrum is only a power- 
law for ^ > ^0, in other words 

P,a{k) = C {kl + ky'" = C {kl + k^p-''^\ (5) 

To a power spectrum of this form corresponds a correlation 
function of fluctuations: 

^(.)=A'(^)"V,/,(^), (6) 

where, rc = lir/kQ, Kjf{x) is the r/-order, modified Bessel func- 
tion of the second kind (also sometimes referred as hyperbolic 
Bessel function), and A = C'2i-''7r-<''+i>/2r?/r[(r/+ l)/2]. 
The Nth-dimensional power spectrum index for a correlation 
function of the form B cx {r / rc)^I^K,-ii2{2iTr / rc) can also be 
generalized to Pnd oc {kl + k^y"!"^, with 

lND = -N-7i. (7) 

This relation is very similar to that for the long wave domi- 
nated case (eq0). Actually for r ^ , Kj^jo can be expanded 
as ~ (r/rc)^l^, thus the 3D correlation function (as opposed to 
the structure function as in the long wave dominated regime) 
goes as a power-law B(r) ~ (r/rc)'' for small separations. No- 
tice that for a shallow spectrum the structure function grows 
rapidly at the smallest scales and then flattens. Similarly to the 
steep spectrum, we produced a short-wave-dominated Gaus- 
sian 3D field, with a prescribed index of 730 = -2.5. In Figure 
|2]we present the expected, and calculated two point statistics. 
Here the critical scale rc is determined by the smallest wave- 
number (feo = 27r/L), in our case it corresponds to the size of 
the computational box {jc = L). We see again a fair agreement 
with the prescribed and calculated spectra. However, Figure 
|2]reveals a significant departure of the calculated correlation 
function with the prediction from equation for large lags. 
The explanation of such difference is that the analytical re- 
lations of spectra (eqs. Q|5l) with structure and correlation 
functions (eqs. ||3]|6l) is exact only in the limit of continuous 
integrals over infinite wave-numbers. The data-sets presented 
in this section are constructed in Fourier space, then trans- 
lated to real space by means of discrete Fourier transforms of 
the form: 

L-l 

U{x) = ^ \Pw(k) I exp [iTTkx/L] (8) 

k=\ 

The sum runs from k= \ ensuring that (m(x)) = 0. In practice, 
we evaluate the Fourier Transforms via FFT and set explicitly 
the k = Q component of the spectrum to zero to guarantee that 
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Fig. 2. — Same as Fig. 1 but for a short-wave dominated spectrum (3D 
spectral index of 730 = -2.5). 



the average of the fluctuating part of the field is null. The re- 
sulting field has a limited range of harmonics available, deter- 
mined basically by the computational grid size (L). We have 
constructed large ID fields, in which we see that the gap be- 
tween the analytical and the computed correlation functions 
gets smaller as we increase resolution. Thus, it is not sur- 
prising that spectra shows a much better correspondence than 
correlation functions. At the same time, structure functions do 
not deal wi th the lowest harmonic s, which introduce largest 
errors (see Monin & Yagloml|l975h . And therefore, they are 
less affected by the lack of lower harmonics as can be con- 
firmed by the fact that they are closer to the analytical predic- 
tion than correlation functions. It is important to always keep 
in mind the issues that can arise from the discrete nature of the 
data. However, we must note that this particular problem lies 
within the generation of the data-sets in frequency space and 
not in the computation of correlation or structure functions via 
spectral methods. We obtain identical results using FFT and 
looping in directly in real space to do the average required. In 
real life, the limitation is likely to be in the opposite direction: 
the finite wave-numbers available would show up as uncer- 
tainty in determining the power spectra while structure and 
correlation functions should be estimated with smaller errors 
(if measured directly in real space). 
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2.2. Structure functions of quantities projected along the 
line of sight 

From spectroscopic observations we can not obtain either 
the density or velocity fields in real space (x, y, z), but we have 
to deal with projections along the line of sight (LOS). Despite 
the fact that our main goal is to extract velocity statistics from 
the centroids of velocity, we will discuss in this section the 
statistics of density integrated along the LOS (column den- 
sity). It will become clear later that the same procedure can 
be applied to obtain velocity statistics . The issue of p rojec- 
tion has been previously discussed in iLazarianl Jl995li) . here 
we briefly state some results that are relevant to this work. In 
Appendix B we exemplify the projection effects of structure 
functions for the particular case of power-law statistics. And 
in Appendix C the power-spectrum of a homogeneous field 
that has been integrated along the LOS. 

In what follows we will assume that the emissivity of our 
medium is proportional to the first power of the density (this 
is true, for instance, in the case of H I). We will consider an 
isothermal medium, and neglect the effects of self-absorption. 
In this case the integrated intensity of the emission (integrated 
along the velocity coordinate) is proportional to the column 
density (see Appendix A for more details): 



/(X)= / ap,(X,v,)dy, 



a p(x) dz, 



(9) 



where a is a constant and p(x) is the mass density. The den- 
sity of emitters /9s(X, v,) can be identified as the column den- 
sity per velocity interval, commonly referred as d^V / dv. To 
distinguish between 2D and 3D vectors, we will use capital 
letters to denote the former and lower case for the latter (i.e. 
X = [x,y], X = [.JCj^jZ]). Our assumption is satisfied for obser- 
vational data where the medium is optically thin, thermalized, 
and with constant excitation conditions. However, for any 
observed map its applicability has to be examined carefully. 
Even for H I widespread self-absorption has been detected (for 
example Jackson et al. 2002; Li & Goldsmith 2003). 

Consider the structure function of the integrated intensity 
described in equation (|9j 



([/(Xi)-/(X2)]') = 




a p(xi) dz- 



a p(x2) dz 



where we have written explicitly the limits of integration, with 
Ztoi being the size of the object (in the LOS direction). Clearly 
Ztot does not necessarily have to coincide with the transverse 
size (in the plane of the sky) of the object under stud y, how- 
ever th at is the case in our data-sets. As described in lLazarianI 
il995h . we can expand the square in equation ilOl combining 



X(x)djc 



and the elementary identity 



x(-*^i)x(^2)dJCidJC2 



(11) 



(a — b)(c—d) = 



to obtain 



1 r 



(a-df + {b-cf-(a-cf-(b-df 



■ a- I I dzidz2 
/o Jo 



([/(X,)-/(X2)]')-"2 

Where dp{r) is the 3D structure function of the density. 
dp(r)={[p(xi)-p(x2)f). 



(12) 



<ip(r)-c/p(r)|x,^X: 
(13) 



(14) 



This definition is general and does not require any particular 
functional form of the 3D structure functions (i.e. power-law 
statistics). The problem of formally inverting equation il3\ to 
obtain the underlying statistics, allowing for anisotropic tur- 
bulence, with an arbitrary spectrum (i.e. not a power-law), has 
been presented in Lazarian ( 1995), but it is somewhat mathe- 
matically involved. For 3D fields with a power-law spectrum, 
homogeneous and isotropic, the structure functions of the in- 
tegrated fields (2D maps) can be simply approximated by two 
power-laws, one at small the other at large separations (see 
Minter 2002, and also Appendix B). For instance, if the den- 
sity has a power-spectrum P3o_p cx k''^", the structure function 
of column density will have the form 



([/(Xi)-/(X2)]2)cx7;^ 



(15) 



where R is the separation in the plane of the sky (R = |R| = 
IX2-X1I), and 

{-73D - 2 for R ^ Ztot (either steep or shallow spectrum), 
-73D - 3 for > Ztot , and 73^ < -3 (steep spectrum), 
for R '3> Ztot, and 73/5 > -3 (shallow spectrum). 

(16) 

In contrast, the power spectrum of a field integrated along 
the LOS corresponds to selecting only the kios = compo- 
nents of the underlying 3D spectra, more precisely only the 
solenoidal part (see Appendix C). Thus, for isotropic and ho- 
mogeneous power-law statistics the 2D power spectrum will 
reflect the 3D spectral index. If for instance the density has a 
power spectrum Pio.p oc k''^", the spectrum of column density 
will scale as 

PiDj^K^'"- (17) 

We computed spectra and structure functions for the Gaus- 
sian cubes used before, integrated along the z direction. And 
presented them in Figure |3| along with the the expected be- 
havior from eas.( ll6t and Hit . For the structure functions 
we show only the asymptotic behavior for small separations 
(R <C Zioi) as the R ^ Ztot scales are unavailable (the maximum 
scale not affected by wrap-around periodicity is L/2 = z,oi/'2)- 

If the LOS size of the object under study is much smaller 
than its size in the plane of the sky (so R ^ Ztot is possible) 
we would see the underlying 3D spectral index of the struc- 
ture functions for large separations (i.e. no projection effect). 
With enough resolution we could use the 2D spectral index 
of the column density map for R ^ Ztot to infer the underly- 
ing 3D index of the density. However, for the resolution used 
here (216^ pixels), the projected structure function for small 
lags is already in the transition between the two asymptotics 
m eq. fT6l . Thus the measured index of the projected map 
is always shallower (smaller) than the actual p w -73D - 2. 
Another subtle, yet interesting point, is that power spectrum 
applied to map of integrated velocity field, such as 



K(X) = / v,(x)dz 



(18) 



will recover only the incompressible (solenoidal) component 
of the field. This could be potentially used to study the role of 
compressibility in turbulence statistics, by combining velocity 
centroids and VGA (LE03). 

3. MODIFIED VELOCITY CENTROIDS REVISITED. 

Velocity centroids have been widely used to rela te their 
statistics with velocity, their conventional form is jMunchI 
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Fig. 3. — Two point statistics for Gaussian fields integrated along the z 
direction. In panel (a) we show the 2D power spectra, panel (b) correspond to 
the second order structure functions. The stars correspond to a steep spectrum 
with a prescribed 3D power spectrum index of 730 = -11/3, the diamonds 
correspond to a shallow spectrum with a prescribed 3D power spectrum index 
of 73D = —2.5. In solid and dotted lines respectively, we plotted for reference 
the expectations for both long-wave or short-wave dominated cases. The 
vertical scales in panel (fc) have been modified arbitrarily for visual puiposes. 



C(X) = 



(19) 



We will refer to this definition as "normalized" centroids. The 
denominator in equation M9\ introduces an extra algebraic 
complication for a direct analytical treatment of the two-point 
statistics. For the sake of simplicity we start considering "un- 
normalized" velocity centroids: 



5(X) = 



a V, p.s(X, V;) dv. 



(20) 



(notice that they have units of density times velocity as op- 
posed to velocity alone). 

Similarly to the expression in equation with emissivity 
proportional to the first power of the density, and no self ab- 
sorption, the structure function of the unnormalized velocity 



([5(Xi)-5(X2)]2) = ^(^a Jv,(xi)p{xi)dz-a J v^ixj) pixj) dz 

(21) 

Presenting the density and velocity as a sum of a mean value 
and a fluctuating part: p = po + p, v, = vo + Vj. Where po = (p), 
vo = (v.), and the fluctuations satisfy (p) = 0, (v-.) = 0. Analo- 
gous to equation ( I13t . the structure function of the unnormal- 
ized centroids can be written as 



{[S(Xi)-S(X2)f)=a' J I dzidzi [D(r)-D(r)|x,=xJ 



with 

Z)(r) = ([v,(xi)p(xi)- v,(x2)p(x2)]') . 
And Z)(r) can be approximated as: 



(22) 
(23) 



D(r) « (v2) dp{r) + dAr) - -dp{r)dAr) + c(r), (24) 

which includes the underlying 3D structure function of the 
LOS velocity 



d„_{r) = {[v,{xi)-v,{x2)f), 



(25) 



and cross-correlations of velocity and density fluctuations "*: 

c(r) = 2 (v,(xi)p(x2))' -4po (p(xi)v,(xi)v,(x2)) . (26) 

Because the derivation of equations ( I22> . ( I24> . and M6\ in- 
volves some tedious algebra, we place it in Appendix D. With 
all of this, the structure function of unnormalized velocity 
centroids can be decomposed as 

([5(Xi)-5(X2)]') =/l(R) + /2(R) + /3(R) + /4(R), (27) 
where 



/2(R) = a2(p2) 

/3(R) = --a2 
2 



dzidz2 [t/p(r)- t/p(r)|x,^x,J ' (28) 
dzidz2 t/„(r)-t/v.,(r)|x,=x2 ' (^9) 
dzidz2 dp(r)d,ir)- dp{r)\^^^^^ d„Xr)\^^^Q 3) 



/4(R) = a2 //dzidz2[c(r)-c(r)|x,=xj 



(31) 



With this new decomposition is more evident the defini- 
tion of the structure function of "modified" velocity centroids 
(MVCs), which is 

M(R) = ( [5(Xi) - S{X2)f -{vf) [/(Xi) - /(X2)]') 

= ([5(Xi)-5(X2)]')-/l(R) 

= /2(R) + /3(R) + /4(R). (32) 

The velocity dispersion (v?) can be obtained directly from ob- 
servations using the second moment of the spectral lines 



2\ _ /^z Ps(X,v.) dv. 



/ p.5(X,v.) dv- 



(33) 



Thus also /I (R), which can be related to the structure function 
of column density as /1(R) = (v^) ([/(Xi)-/(X2)]2). 

Note that equation 1261 is somewhat different from LE03, where there 
was a misprint; which has no effect on the results presented. 
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Similarly, the power spectrum of centroids can be decom- 
posed as (details in Appendix F): 

^2D,s(K) = {vl)(az,o,f6(K) + v^P2d,/(K) + a^plP2D.v-SK) 
+J^{B3(R)} + J^{B4(R)}. (34) 

The term {p'^){v^)(aztot)^S(R) has no effect in the slope 
of the power spectrum because it only has power at K = 0. 
P2D./(K), and P2d.V:(K.) are the spectra of column density, and 
integrated velocity respectively. They can be used to obtain 
the 3D spectral index of density (or velocity) as shown in 
Appendix C. J-{B3(R)} is the Fourier transform of B3(R), 
a cross-term term analogous to /3(R), but in terms of cor- 
relation functions. Similarly, JF{Z?4(R)} include the same 
density-velocity cross-correlations as /4(R). 

The power spectrum of MVCs can be obtained by subtract- 
ing (vJ)P2D./(K) from the spectrum of centroids. We derive in 
Appendix G a criterion for MVCs to trace the statistics of ve- 
locity better than unnormalized centroids. It was found that, 
with very little dependence on the spectral index, MVCs are 
advantageous compared to unnormalized centroids at small 
lags. This result is general and we tested it using analytical 
expressions for the structure functions. For simplicity we con- 
sidered only the two physically motivated cases, steep density 
with steep velocity, and shallow density with steep velocity. 
The latter was not included explicitly in Appendix G because 
we obtain almost identical results in both cases. 

If v'o = and T {B3(R)} + T {B4(R)} can be neglected the 
spectrum of unnormalized centroids will trace the solenoidal 
component of the underlying velocity spectrum (FffN[K,Q], 
see Appendix C). But if the turbulent velocity field is 
mostly solenoidal, as supported by numerical simulations 
jMatthaeus et al. 1996; Porter, Woodward & Pouquet 1998; 
iCho & Lazarianll2003l) . the power-spectrum is uniquely de- 
fined assuming isotropy (E{k) = J PND(^)dk w 4Trk^FNN[k]). 
In the same way if /1(R) can be eliminated (either for being 
small compared to the structure function of centroids or by 
subtraction -MVCs-) and if /2(R) > /3(R) + /4(R), the struc- 
ture function of the remaining map will trace the structure 
function of a map of integrated turbulent velocity. And we 
can in principle recover the underlying 3D velocity statistics 
(see Appendix B). With this background (and disregarding the 
cross-terms /3[R], /4[R]) we arrived in LE03 to a criterion for 
safe use of (unnormalized) velocity centroids: ;/ 

{[S(Xi)-S(X2)f) » (v2) ([/(Xi)-/(X2)]2) , (35) 

then the structure function of velocity centroids will mostly 
trace the turbulent velocity statistics, otherwise the density 
fluctuations are important and will be reflected in the centroid 
measures. When the structure function of velocity centroids 
is shallower or at least not much steeper than that of the col- 
umn density, which can be verified by the power spectrum, or 
computing the structure function directly in 3D for a few of 
values for the lag. Then the criterion proposed in LE03 can 
be simplified to use only the variances of the two maps (and 
the velocity dispersion): 

{S'-)»{^^){P). (36) 

If any of these two criteria is violated, one could in principle 
subtract the contribution of density and the MVCs would trace 
velocity structure function, provided that we could neglect the 
cross-terms. 

The contribution of velocity-density cross-correlations 
(c[r]) have been studied earlier. For VCA it has been shown to 



be marginal JLazarian et al l2()()!llEsauivel et alJ2003h . How- 
ever, a more detailed discussion of their effect in the context 
of MVCs is necessary, and is provided below. 

/3(R) and JF{B3(R)} are in some sense "cross-terms", 
/4(R) and T {BA{R)}) are related to correlations between 
density and velocity. We expect both pairs to grow as we 
increase the "interrelation" between density and velocity. We 
will refer to /3[R] (or JF{B3(R)}) simply as "cross-term", and 
to /4[R] (or JF{Z?4(R)}) as "cross-correlations" of density- 
velocity. The latter should be zero for uncorrelated data. At 
the same time the cross-term can be studied analytically for 
power-law statistics as presented in Appendix E, and will 
not be zero, even in the case of uncorrelated velocity and 
density fields. Before computing them directly, in order to 
get a feeling of how important the cross-term could become 
one can consider structure functions. First of all, note that 
([5(Xi)-5(X2)]^) is positive defined, and so are /1(R) and 
/2(R). The remaining terms can be negative, in which case 
they must be smaller than the sum of /1(R) and /2(R). Let us 
focus for the moment on the contribution of /3(R) and dis- 
regard cross-correlations between density and velocity. Its 
magnitude is maximal at large scales, and so are /1(R) and 
/2(R). At such scales |/3(R)| is on the order of /1(R) and 
/2(R). However, in /2(R) the structure function of velocity is 
weighted by {p^) = Pq+ (p^) instead of only (p^), enhancing 
the velocity statistics compared to the cross-term. The im- 
portance of the cross-term at the small scales (in which we 
are most interested) will depend on details such as how steep 
the underlying structure functions are (see Appendix E), and 
the zero levels of density and velocity (see Osse nkopf et al.l 
2005). This is easy to understand for a particular case of 
power-law statistics of the form dp{r) cx r", and c/,.„(r) cx r™, 
with (;«, n > 0, i.e. both fields steep). Here the cross- 
term scales as oc r'""*"", steeper than both velocity and den- 
sity. If at large scales /1(R) and /2(R) are on the order 
of |/3(R)|, provided that the latter falls more rapidly toward 
small scales, its contribution will be smaller than both ve- 
locity and density structure functions, at those scales. But 
if the density or the LOS velocity (or both), have a shal- 
low spectrum, the cross-term can be larger than /1(R) or 
/2(R), and can affect significantly the statistics of centroids. 
Measured spectral indices of density in the literature, range 
from 73/5 ^ -2.5 to 730 ~ -4.0, which include both shal- 
low and steep. This is true for obse rvations in different en- 
vironments in the ISM (for instance, Deshoande et al. 200(1 
Bensch, Stutzki & Ossenkocf 2001; Stanimirovic & Lazaria^ 
2001; Ossenkocf & Mac Low 2002), as well as for numeri- 
cal simulations (see C ho & Lazarian 2002; Bruiit & Mac Lowl 
2004; Beresnvak et alj l2005l) . The velocity spectral index is 
less known from observations, but has been measured to be 
only in the steep regime (for example using VCA, e.g. in 
Stanimirovic & Lazarian 2001), also in agreement with simu- 
lations. From the theoretical standpoint, at small scales when 
self-gravity is important we might expect clumping that result 
on enhanced small scale structure (yielding a shallow spec- 
trum). On the other hand there are no clear physical grounds 
to our knowledge that will produce a small scale dominated 
(shallow) velocity field. However, even in the simple case 
of steep density and steep velocity spectra it is not clear be- 
forehand how important density-velocity cross-correlations 
(74 [R]) could be. Later, we will analyze the contribution 
of the cross-terms and density-velocity cross-correlations in 
more detail, including spectra. 
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TABLE 1 

Parameters of the four runs used. 



Model 


Bo 


(Pgas) 


/3 


rms V 


Ms 


Ma 




A 


1.0 


2.0 


4 


<0.7 


~0.5 


~0.7 


93.9 


B 


1.0 


0.1 


0.2 


^0.7 


~2.5 


~0.7 


4.6 


C 


0.1 


0.1 


0.2 


<0.7 


~2.5 


~ 8 


4.2 


D 


1.0 


0.01 


0.02 


^0.7 


~7 


~0.7 


2.2 



4. TESTING VELOCITY CENTROIDS NUMERICALLY 

In LE03 we performed some preliminary tests of the modi- 
fied velocity centroids using numerical simulations and com- 
pared the power-spectrum with that of velocity field, normal- 
ized (equation 1191 ). and unnormalized centroids (equation 
1201 ). In this section we provide a more detailed test to inves- 
tigate under what conditions velocity centroids can be used to 
recover the velocity statistics. 

4. 1 . The data 

We took compressible MHD data-cubes from the numeri- 
cal simulations of Cho & Lazarian (2003). This data-cubes 
correspond to fully-developed (driven) turbulence. The tur- 
bulence is driven in Fourier Space (solenoidally) at wave- 
numbers 2 < (k^rivingL / [2Tr]) < 3.4. The data-cubes have a 
resolution of 216^' pixels. We use four sets of simulations, the 
parameters for each run are summarized in Tabled The 
models include various values of the plasma (3 (ratio of gas 
to magnetic pressures), sonic Mach numbers j, and Alfven 
Mach number M^- All of these parameters can be found in 
the ISM under different situations. For mo re details about 
the simulations we refer the reader to lCho &'Lazarianl2003.) . 
The outcome of the simulations are density and velocity data- 
cubes that we use to compute the centroids. We will refer to 
this data-sets as "original", and the existing correlations be- 
tween density and velocity (consistent with MHD evolution) 
are left intact. The numerical simulations have a limited iner- 
tial range. We do not have power-law statistics (i.e. inertial 
range) at the largest scales (smallest wave-numbers) due to 
the driving of the turbulence. And neither we have power- 
law statistics at the smallest scales (largest wave-numbers), 
because of numerical dissipation. Thus, it is very difficult 
to estimate spectral indices because the measured log-log- 
slope is quite sensitive to the range in wave-numbers (or lags) 
used. This poses a problem of obtaining quantitative results. 
For that reason, we created another data-set by modifying the 
original fields to have strict power-law spectra, following the 
procedure in Lazarian et al. (2001). The procedure consists 
in modifying the amplitude of the Fourier components of the 
data so they follow a power-law, while keeping the phases in- 
tact. This way we preserve most of the spatial information. By 
keeping the phases we also minimize the effect of the modifi- 
cation to the density-velocity correlations. In addition, these 
new fields have the same mean value than the original data- 
sets, and the magnitude of their power-spectra (vertical offset) 
was fixed to match the original variances as well. We will re- 
fer to this data-sets as "reformed". 

4.2. Results 

Column density and centroids of velocity are two- 
dimensional maps, and therefore it is not computationally re- 
strictive to obtain their correlation or structure function di- 
rectly in real space. Power-spectrum is often computationally 



TABLE 2 

Three-dimensional spectral indices. The values 
IN parentheses correspond to the original 

SIMULATIONS, AND IN BOLD FACE TO THE reformed 
DATA-SETS. 



Model 


Density 


LOS Velocity 


-73D "1 


-73D 


A 


(3.5) 3.5 (0.4) 0.6 


(3.8) 3.8 (0.5) 0.8 


B 


(3.3)3.3 (0.3)0.4 


(3.6) 3.6 (0.5) 0.6 


C 


(3.1)3.1 (0.1)0.3 


(4.0) 4.0 (0.8) 0.9 


D 


(2.6)2.6 ■■■ * 


(3.8) 3.8 (0.6) 0.8 



The measured power-spectrum index in this case conesponds to a 
shallow spectrum. Thus, the conelation function is expected to follow 
a power-law, not the structure function. 



cheaper because EFT can be used. However, inherent difficul- 
ties of applying Fourier analysis to real data, for instance the 
lack of periodic boundary conditions and instrumentational 
response make power spectrum often unreliable. To allevi- 
ate this problem more elaborate techniques like wavelet trans- 
forms have been proposed (Zielinskv & Stutzki 1999). More- 
over, if the observed maps are not naturally arranged in a 
Cartesian grid one would need to smear the data onto that 
kind of grid to use FFT. Which is not necessary for struc- 
ture or correlation functions if computed directly averaging 
in configuration space. In despite of this, because the simu- 
lations we used are in a Cartesian grid and indeed have pe- 
riodic boundary conditions, we can compute spectra, corre- 
lation, and structure functions using FFT (see §2) with as 
good accuracy than doing the average in real space. The rela- 
tion between the structure function, correlation function, and 
power spectrum of unnormalized centroids can be found in 
Appendix F. 

4.2.1. 3D statistics 

Before we study the 2D maps and try to extract from those 
the underlying 3D statistics we will start by computing the 
three-dimensional statistics. This is shown in Figure |3 It is 
noticeable that only for the case where the turbulence is sub- 
sonic {M.S ~ 0.5) the level of the velocity fluctuations is larger 
than that of the density. In the figure is also evident the lim- 
ited inertial range in the original simulations (i.e. not perfect 
power-law statistics). Spectral indices (log-log slopes), both 
for power-spectra and structure functions from Figure |3 are 
given in Table |2] For the original data-sets the indices 

for power-spectra were obtained in a range of wave-numbers 
kL/{2Tr) ~ [5- 15] (between the scale of injection and the 
scales at which dissipation is dominant). The structure func- 
tions spectral indices were calculated with the corresponding 
values for spatial separations, r/L ^ [1/15- 1/5]. The re- 
formed data-sets were constructed using the power-spectra in- 
dices estimated for the original MHD simulations. We can see 
in Fig. |4]the idealized power-law spectra of the reformed data- 
sets. At the same time, although structure functions do not 
have such perfect power-law behavior (see discussion at the 
end of §§2.1), they show improvement compared to the origi- 
nal simulations. The range of scales used to measure the spec- 
tral indices for the reformed data-sets is kL/ilii) ^ [3.5-100] 
or r/L ~ [1/100- 1/3.5], much wider than for the original 
sets. Notice that power-spectra for the density becomes shal- 
lower with the increase of the Mach number At the same time 
the spectral index of velocity is always steep. 
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Fig. 4. — Underlying 3D statistics of the MHD simulations. At the top the power spectra, and on the bottom the corresponding structure functions. The four 
runs are ordered in ascendant Mach number from left to right. 



4.2.2. Statistics of projected quantities (2D) 

A natural way to study how the velocity centroids trace the 
statistics of velocity is to compare their two point statistics 
to those of an integrated velocity map ("equation II 8>. This 
map can be used to obtain the velocity spectral index in the 
same way column density can be used to obtain that of den- 
sity. Therefore it is a direct measure of the underlying velocity 
statistics (see Appendices B and C). However, it is not ob- 
servable, while velocity centroids are. We computed power- 
spectra and structure functions of 2D maps of the various 
centroids (normalized, unnormalized, and "modified"), inte- 
grated velocity, and integrated density. The results for power- 
spectra are shown in Figures |5j and |6j for the original, and 
modified data-sets respectively. Similarly, Figures and |8l 
show the results for structure functions. The most notice- 
able difference in the Figures is of course the larger inertial 
range of the reformed data-set. We also present in Table |3l 
a summary of spectral indices (log-log slope) measured over 
5 < KL/(2tt) < 25 (or 1/25 <R/L< I /5) for the original data 
sets, and kL/iln) ~ [3.5-100] (or r/L - [1/100- 1/3.5]) for 
the reformed sets. 

Comparing the spectral indices derived in 3D (Table|2} with 
those of column density and integrated velocity in Table |3j 
one can notice a better correspondence for the reformed data- 
sets. This is true for the power-spectra index 73/3; as well as 
for m, and /i for structure functions (related by equation ll6l ). 
Directly from Figures (SHS] one can see that only for the case 
of subsonic turbulence (model A, Ms ~ 0.5) the spectrum of 
centroids clearly scales with that of integrated velocity. In 
this case the power-spectra of all the variations of centroids 




K L/(2n 



K L/{27i ) 



Fig. 5. — Power spectra of the integrated density (crosses), integrated ve- 
locity (stars), unnormalized, normalized, and modified centroids (solid, dot- 
ted, and dashed lines respectively) for the original set of simulations. We 
multiplied the spectrum of velocity fluctuations by pj, and that of normalized 
centroids by (/- ) . We show the spectrum of integrated density for reference 
only, to be in the same units as the other quantities in the figure it should be 
multiplied by Vq, but since vg ~ the scaling is omitted here. 
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TABLE 3 

Spectral indices (2D) of quantities integrated along the LOS. The values in parentheses correspond to the 

original SIMULATIONS, AND IN BOLD FACE TO THE reformed DATA-SETS. 





Integrated 


Integrated 


Unnomialized 


Normalized 


Modified 


Model 


Density 


Velocity 


Centroids 


Centroids 


Centi'oids 




-73D 


-73D 


-73D 


-73D 


-73D 


A 


(4.0) 3.5 (0.5) 1.4 


(3.6) 3.8 (0.9) 1.5 


(3.6) 3.8 (0.8) 1.5 


(3.5) 3.8 (0.8) 1.5 


(3.6) 3.7 (0.8) 1.5 


B 


(3.8) 3.3 (0.4) 1.2 


(3.9) 3.6 (I.O) 1.4 


(3.5) 3.3 (0.8) 1.1 


(3.5) 3.2 (0.9) 1.1 


(3.8) 3.3 (1.0) 1.0 


C 


(3.3) 3.1 (0.6) 1.1 


(4.5)4.0 (1.3)1.6 


(3.8)3.4 (1.1)1.3 


(3.9)3.5 (1.2)1.3 


(3.7)3.4 (1.2)1.5 


D 


(2.8) 2.7 (0.2) 0.7 


(4.7) 3.8 (0.8) 1.5 


(3.4) 2.8 (0.5) 0.9 


(3.8) 2.9 (0.7) 1.0 


(3.5) 2.9 (0.8) 1.1 



Since this index is not measured at scales con'esponding to ^ Ziot, but rather in the transition between the two asymptotic regimes in eq. lI6l . this are only 
lower hmits on the actual /x for small lags. 




Fig. 6. — Same as Fig.|5] for the reformed data-sets. 



recover the spectral index of velocity, within 10% error for 
the original simulations, and < 3% for the reformed data sets. 
For the cases of mildly supersonic turbulence (models B, and 
C, M.S ^ 2.5) is not clear neither from the Figures, nor from 
the measured slopes. While for the strongly supersonic case 
(model D, Ai^ ~ 7) it is obvious that velocity centroids fail to 
recover the velocity scaling. 

Due to the finite width effects discussed at the end of §§2.2, 
it is more difficult to determine quantitatively the spectral in- 
dex from centroid maps using structure functions. According 
to equation M6\ . one should restrict to measure the spectral 
index for lags either much smaller, or much larger than the 
LOS extent of the object under study. The latter is not feasi- 
ble with our simulations because the maximum lag available, 
unaffected by wrap-around periodicity, is L/2 = Ztot/'^- The 
other case (R <^ Ztot) is not strictly possible with the resolu- 
tion used here. For the MHD data-cubes we have to avoid the 
smallest scales because they are not within the inertial range 
(i.e. they are already dominated by dissipation). Actually, 
the lags used to measure fi in the original data-sets are well 
in the transition between the two asymptotic power-laws of 
equation M6\ . Thus, if we obtain the 3D spectral indices as 



Fig. 7. — Structure functions of integrated density (crosses), integrated 
velocity (stars), unnormalized, normalized, and modified centroids (solid, 
dotted, and dashed lines respectively) for the original set of simulations. . 
To allow for a direct comparison we scale the quantities to be all in units of 
[v^p^] (i. e. we multiplied the structure function of the integrated velocity by 
(p"), the integrated density by (i'-), and the normalized centroids by (/')). 

73£) ~ -/i-2 (taking /i from Table|3}, we would only get upper 
limits of the actual 730. Keeping in mind that our definition 
of the index includes the minus sign, this means that the real 
index is going to be more negative than the inferred 73/5 from 
structure functions. A lower limit on the spectral index can be 
obtained as 730 ~ -/x-3, unfortunately this provides a rather 
wide range of possible indices and it is not useful for practical 
purposes (for instance distinguishing between two models of 
turbulence). For the reformed data-sets the situation is better, 
because we can use the smallest scales. This can be verified 
from the better agreement of 730 ~ -/i-2 with 73/3 measured 
directly from power-spectra. However, R <C Ztot is still not 
strictly fulfilled, and thus we get only lower limits of /i. We 
recognize this projection effect as an important drawback for 
the structure functions compared with spectra. Power spectra 
do not suffer so strongly of finite width projection effects, and 
could be used to obtain more accurately the spectral index of 
the underlying 3D index fie ld from integrated quantities (see 
also lOssenkopf et'al]l2005l) . This is true for synthetic data. 
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Fig. 8. — Same as FigQ for the reformed data-sets. 

but for real observations power spectrum is not as reliable (see 
iBensch et al. 2001 ). At the same time, the problem of recover- 
ing quantitatively the spectral index from structure functions 
might be overcome for real observations with a larger iner- 
tial range. With enough spatial resolution one can have eas- 
ily over two decades of inertial range in the plane of the sky 
(below the size of a given object along the LOS direction). 
In the spirit of keeping the description in terms of structure 
functions, which might be better suited for some type of real 
observations. And, using the limited range we have, it is still 
possible to do a comparative analysis between the spectral in- 
dex of structure functions of centroids to that of integrated 
velocity or density. 

From Figures|5HElwe see that the relative importance of the 
density and velocity statistics on the centroids changes with 
the Mach number of turbulence. For the low M.s model (pan- 
els [a] in Figs. [SHE} the level of the density fluctuations is very 
small compared to the velocity fluctuations (weighted by (vj) 
and {p^) respectively) and all the centroids trace remarkably 
well the velocity statistics. In this case the simplified criterion 
{S^)/{{v\){P)) > 1 is well satisfied: {S^) / ({vfi {1^)) > 90 for 
the original data-set, and > 30 for the reformed one. This re- 
sult also agrees with the notion that velocity centroids trace 
the statistics of velocity in the case of subsonic turbulence 
(where density fluctuations are expected to be negligible, as 
the turbulence is almost incompressible). For the rest of the 
models density has an increasing impact which is reflected in 
the centroids, and our criteria is either not entirely met, or vi- 
olated. We can also see for supersonic turbulence the slope 
of the structure function of unnormalized centroids is almost 
always steeper than that of column density. This means that 
one need to use the full criterion as proposed in LE03 to judge 
whether centroids will trace velocity or not. We will see how- 
ever that our simplification of the criterion seems to discrimi- 
nate when centroids trace the statistics of integrated velocity, 
at least for the data sets employed here. For all the supersonic 
cases the measured power-spectrum index from the the differ- 



ent centroids fail to give the index of velocity. And in general, 
their centroids index is found to lie between those of veloc- 
ity and density, with more scatter for the original simulations 
compared to the reformed fields. For the original data-sets of 
models B, C, and D the ratio (S^) / {{vD (i^)) is w 5,4,and3 
respectively. For the corresponding reformed data-cubes our 
simplified criterion gives {S^) / {{vT) {r)) k, 2, 4, and 3. In all 
of these cases {S^) /((v?) {r)) ^ 1 is not strictly true, and in- 
deed there is an important contribution of density fluctuations. 
This is consistent with the results in LE03 where we see evi- 
dence of a density dominated regime at small scales (rather 
density contaminated regime as centroids do not show the 
same index as density either). Our results agree with those 
presented in lBrunt & Mac Lo^ 1200 4), where they found that 
centroids do not provide good velocity representation for the 
supersonic turbulence they studied (Al, > 1.9). 



4.2.3. Cross-terms and density-velocity cross-correlations 

The cross-terms, /3(R), JF{B3(R)}, as well as the cross- 
correlations in /4(R), and JF{Z?4(R)} have contributions of 
velocity and density that we cannot disentangle entirely from 
observables. Furthermore, they cannot be expressed in terms 
of integrated quantities and they have to be computed from 
3D statistics. We used Fourier transforms to obtain the struc- 
ture, and correlation functions in 3D needed to produce in- 
dependently all the terms in equations Mil , and ( I34> . These 
3D statistics were integrated numerically to get JF{B3(R)}, 
/3(R), J^{B4(R)}, and /4(R). To check the accuracy of the 
3D quantities obtained and our numerical integration to 2D 
we verified with the cases in which the statistics can be also 
obtained directly in 2D (e.g. /1[R], /2[R]), finding always a 
good agreement. The results of the decomposition in terms 
of power-spectra (equation 1341 . after K average) are shown 
in Figure|9lfor the original simulations, and in Fi sure [Tol for 
the reformed data-sets. The decomposition for R averaged 
structure-functions is presented in in Figures ^2 and^J for 
the original, and modified data-sets respectively. Some- 
thing to notice is that, because structure functions span over 
fewer decades in the vertical axis, the separation of all the 
terms is generally more clear than for power-spectra. It is 
also to be noticed, that for the subsonic model (A), the spec- 
trum and structure function of centroids are almost unaffected 
by density fluctuations, cross-terms, or density velocity cross- 
correlations. Cross-correlations are found to be larger for su- 
personic turbulence compared to the subsonic case, but from 
this limited data set we can not conclude that they scale in 
some particular way with the sonic Mach number. It is also 
significant that the cross-term increases relative to the other 
terms with Ms- For spectra, in all the supersonic models 
(B, C , and D) the cross-term JF{B3(R)} is dominant. This 
would mean that the log-log slope measured from spectra in 
all of this cases is not a direct reflection of the velocity spec- 
tral index. At the same time we observe that the magnitude 
of density-velocity cross-correlations can only be entirely ig- 
nored for model A (subsonic turbulence). For structure func- 
tions the cross-term is comparable or larger in magnitude than 
the contribution of column density. At the same time it gets 
closer, but does not become larger than /2(R). In fact, it was 
found to be always steeper than the velocity term, therefore its 
contribution at the small scales could be neglected. Velocity- 
density cross-correlations in /4(R) for all the cases we con- 
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Fig. 9. — De com position of the spectra of unnormalized velocity centroids 
as in equation 1341 (with vq ~ 0), for the original data-sets. The solid line 
is the spectrum of the 2D map of velocity centroids. The dotted Hne (for 
reference only) is the power spectrum of column density, the dashed is the 
spectrum of integrated velocity, both obtained from integrated (2D) maps. 
The crosses represent the cross-term J?-'{i}3(R)}, and the "X" correspond to 
density-velocity cross-correlations JF{B4(R)}. This last two were computed 
from 3D statistics, then integrated along the LOS. The diamonds show the 
sum of all the terms in the decomposition to compare with the solid line. 



Fig. 1 1 . — Dec omp osition of the structure function of unnormalized veloc- 
ity centroids (eg. 1271 ). for the original data-sets. The solid line is the struc- 
ture function of the 2D map of velocity centroids. The dotted line is /1(R), 
the dashed /2(R), both obtained from 2D maps. The cross-term |/3(R)| as 
crosses, while the density-velocity cross-correlations (|/4[R]|) are the "X" . 
This last two were computed from 3D statistics, then integrated along the 
LOS. The diamonds show the sum of all the terms in the decomposition to 
compare with the solid line. 
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Fig. 12. — Same as ll II for the reformed data-sets. 
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Fig. 10. — Same as Fig.|9| for the reformed data-sets. 



sidered here'' never had an important impact on the statistics 
of centroids. Anyway, /3(R) + /4(R) as a whole (especially 
at the smallest separations, which we are most interested in) 



^ The situations when velocity-density correlations are dominant are dis- 
cussed in our next paper 
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has been found to be smaller than the integrated velocity term. 
Suggesting that there could be cases where /3(R)+/4(R) -the 
cross-term and cross-correlations- can be neglected. How- 
ever, since the indices from power-spectra of centroids failed 
to give the velocity spectral index for models B, C, and D, we 
conclude that the retrieval of spectral indices from velocity 
centroids over the entire inertial range should be restricted to 
very low sonic Mach numbers (Ais < 2.5). Furthermore, the 
relative importance of the cross-terms grow together with the 
column density term as the strength of turbulence increases. 
For Ms ^ 7 certainly the contribution of the column density 
is considerable. In this case MVCs could be used to remove 
the contribution from the column density term in the structure 
function of centroids. However at such high Mach number the 
cross-terms cannot be neglected either. 

At a first glance it might seem surprising that the velocity 
centroids could be dominated by velocity even when the 3D 
structure functions of density had a zero point (offset) larger 
than that of velocity. This effect arises primarily from the 
fact that the density is positive defined, and necessarily has 
a non zero mean whereas the velocity field can have a zero 
mean. This results in the factors that multiply the density 
and velocity structure functions in equation i24\ . The factor 
(vj) = Vq + (v'^) multiplies the density structure function, and 
we can minimize the undesired density contribution by shift- 
ing our velocity axis in such a way that vq = 0. For a more 
det ailed discussion a bout the velocity and density zero levels 
see lOssenkopf et al] ( 12005 ). 

5. VELOCITY CENTROIDS AND ANISOTROPY 
STUDIES IN SUPERSONIC TURBULENCE 

We have seen how density fluctuations in supersonic tur- 
bulence affect our ability to determine the velocity spectral 
index from observations. But is there something else ve- 
locity centroids could be used for, even in supersonic turbu- 
lence? The presence of a magnetic field introduces a pref- 
erential direction of motion, thus breaking the isotropy in 
the turbulent cascade. In a turbulent magnetized plasma, ed- 
dies become elongated along the direction of the local mag- 
netic field. Velocity statistics have been suggested to study 
this anisotrop v ( Lazarian et al. 2002; Esauivel et al. 2003 
IVestuto et a l. 2003)^. For instance, iso-contours of two point 
statistics of velocity centroids instead of being circular, as in 
the isotropic case, are ellipses with symmetry axis that reveals 
the direction of the mean magnetic field. We present in Figure 
E]contours of equal correlation of velocity centroids (unnor- 
malized) from our simulations. The magnitude of the mag- 
netic field determines how much anisotropy will be present. 
We can see from Figure that the anisotropy is very clear 
for models A, B, and D, regardless of the large differences in 
sonic Mach number and plasma f3. The only case in which the 
anisotropy is not evident (model C) is our only super- Alfvenic 
simulation. The concept of super- Alfvenic tur bulence is ad- 
vocated, for instance, bv lPadoan et alJ (l2004bh for molecular 
clouds. 

Another way to visualize the anisotropy is to plot the cor- 
relation functions of centroids in the parallel or perpendicu- 
lar directions relative to Z?o, this is shown in Figure O In 
our simulations this corresponds to plot the value at the in- 
tercepts in Figure in observations it would mean to plot 

' Both velocity centroids and spectral congelation functions were demon- 
strated to trace the magnetic field in Lazarian et al. 1 2002 ), channel maps 
were stu died for the same pur pose in iEsouivel et aL 1.2003,1 . and velocity cen- 
troids in lVestuto et alj 120031) 
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Fig. 13. — Anisotropy in the congelation functions: contours of equal cor- 
relation for the simulations (solid lines). For reference we show isotropic 
contours as dotted hnes. The sonic and Alfven Mach numbers (Al,,, Ma 
respectively) are indicated in the title of the plots. The anisotropy reveals the 
direction of the magnetic field for all the sub-Alfvenic cases, regardless of 
the sonic Mach number 
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Fig. 14. — Conelation functions taken in directions parallel and perpen- 
dicular to the mean magnetic field. The anisotropy shows in the different 
scale-lengths for the distinct directions. It is noticeable the little difference 
in panels (a), (b), and (d). Which conespond to the same ratio B/Bq, but 
very different sonic Mach numbers. Panel (c) conesponds to super-Alfvenic 
(B > Bo) turbulence, and the anisotropy is not evident in the centroid maps. 

the correlation function only along the major or minor axis 
of the contours of equal correlation, the anisotropy is evi- 
dent as a different scale-length (correlation-length) for corre- 
lations in the parallel or perpendicular to the mean magnetic 
field. For subsonic turbulence the degree of anisotropy re- 
flects the ratio (B / Bq)-. However for supersonic turbulence as 
the contributions of density get important, and as density at 
high Mach number gets more isotropic (see iBeresnvak et alJ 
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l2005h . one could expect the anisotropy to decrease while the 
ratio B/Bq stays the same. We found very little evidence 
of this decrease, as can be verified from models A, B, and 
D in Figure [21 all the three runs have the same B/Bo, and 
a large contrast in sonic Mach number. In model C even 
though the magnetic pressure is larger than the gas pressure 
{j3 = 0.2), it corresponds to a weak mean field B > Bo- In 
this case the magnetic field has very chaotic structure at large 
scales (see Cho & Lazarian 2003). As the fluctuations at large 
scale determine the anisotropy of the projected data, the MHD 
anisotropy is erased after LOS averaging. 

Whether the ISM is sub-Alfvenic is still up for debate. The 
fact that centroids are only anisotropic for sub-Alfvenic tur- 
bulence gives us an opportunity to study the conditions un- 
der which that is the case. It is certainly encouraging, that 
for sub-Alfvenic turbulence the degree of anisotropy does not 
seem to be strongly affected by density fluctuations. How- 
ever, we should add a word of warning in trying to determine 
the ratio B/Bq from velocity centroids in supersonic turbu- 
lence. The fact that centroids do not represent velocity for 
A4s > 2.5 calls for some caution as to what extent the ratio 
B/Bq can be obtained from observations. Nonetheless, the 
result on the direction of the mean field is rather robust, in- 
cluding for highly supersonic turbulence, provided a relatively 
large ordered component Bo- Anisotropics of turbulence mea- 
sured by different techniques, including centroids, provide a 
promising way of measuring the direction of the perpendicu- 
lar to the LOS component of the magnetic field. For practical 
purposes the technique can be tested using polarized radiation 
arising from aligned dust grains (see Lazarian 2003 for a dis- 
cussion of grain alignment), or CO polarization arising from 
Goldreich-Kylatis effect (see lLai. G irart & Crutcher 2003). 

For any particular observational data one should bear in 
mind that only the plane-of-sky components of mean and fluc- 
tuating magnetic field are available. Therefore, for instance, 
a cloud with sub-Alfvenic turbulence but with a mean mag- 
netic field directed along the line of sight would look like 
a cloud with super-Alfvenic turbulence. The distinction be- 
tween the two cases may be done, however, using an ensem- 
ble of clouds. It is unlikely to have mean magnetic field al- 
ways directed towards the observer. If, indeed, the turbulence 
is typically sub-Alfvenic, the anisotropy should show up in 
centroid maps. 

6. DISCUSSION 

We used MHD data-cubes to produce centroid maps, and 
explored the limitations of velocity centroids for studying in- 
terstellar turbulence. 

We investigated numerically the analytical predictions in 
our previous study of velocity centroids (LE03). In that work, 
we decomposed the structure function of velocity centroids 
into three main contributions, namely integrated (or column) 
density, integrated velocity, and "cross-terms". By calculating 
separately all the terms in the structure function of unnormal- 
ized centroids we showed that the decomposition works well. 

We found two important restrictions on the procedure accu- 
racy, both related to the finite resolution in the numerical sim- 
ulations. First, structure functions of integrated quantities can 
in principle be used to retrieve the underlying 3D spectral in- 
dex if calculated for lags either much greater or much smaller 
than the LOS extent of the object. We showed however, that 
with the resolution used here (216^') structure function of in- 
tegrated quantities can only give an upper limit on the actual 



spectral index. Second, the limited inertial range in our sim- 
ulations present an additional constraint on the accuracy of 
estimating the spectral index, mainly because the measured 
log-log slope is very sensitive the range of scales used. 

In regards to the finite width projection effects one can po- 
tentially still do a comparative study of the centroid statistics 
with the statistics of velocity and density. That is, we can 
compare the spectral index of the centroid maps to that of col- 
umn density or integrated density, and see to which one is 
closer. However, for real observational data with the possibil- 
ity to sample lags R <C Zwi, obtaining the 3D spectral index 
directly from structure functions is not a problem. 

Power-spectrum is an alternative to obtain the scaling prop- 
erties of the turbulent velocity field that does not suffer of such 
projection effects (i.e. the spectral index can be measured at 
any wave-numbers within the inertial range). In order to use 
power-spectra too, we constructed a power-spectra decompo- 
sition analogous to that in LE03. And, by calculating each 
term separately, we tested successfully the validity of this de- 
composition. 

We found that the scaling properties of the underlying ve- 
locity field (i.e. spectral index) can be reliably retrieved for 
subsonic turbulence. However, the contamination from cross- 
terms clearly showed up in all the supersonic cases for spectra, 
while for structure functions it was only evident for A^j ^ 7. 
To alleviate the other problem (not enough inertial range) we 
introduced reformed versions of the original simulations. The 
power spectra of new data-set are strict power-laws, with al- 
most the same spatial structure, and density-velocity cross- 
correlations. The spectral indices obtained with these re- 
formed data sets are in better agreement with the analytical 
predictions. 

We also tested successfully a criterion for which velocity 
centroids give trustworthy information (equation 1351 ). that 
can be obtained entirely from observations. For practical pur- 
poses we suggest as the first approach an approximation of the 
criterion in LE03 in terms only of the variances of the maps 
of column density and velocity centroids instead of computing 
all the structure function (equation 1361 '). If this ratio is less 
than unity velocity centroids are not trustworthy. When the 
ratio in ea.(l36> is small but larger than unity it is recommend- 
able to use the full criterion, as proposed in LE03. In our only 
case of subsonic turbulence (model A) centroids traced the 
integrated velocity structure function extremely well. And, 
our simplified criteria was fulfilled: {S^)/{{vl){P))>90, 30, 
for the original, and reformed data-set respectively. For the 
rest of the models (B, C, D) the criteria was not satisfied 
({S-) /{{vfi (/-)) > [5-2]), and the centroids failed to trace the 
statistics of the velocity field. In this models one can see a bet- 
ter correspondence between centroids and integrated velocity 
for the Ais ^2.5 runs compared to the highly supersonic case 
Alj ^ 7. Thus, one might hope centroids to be able to trace 
the spectral index of velocity for supersonic turbulence, but 
only for Ai^ < 2.5. We must recognize that this study was 
done using a limited data set. A more complete exploration of 
parameter space is desirable to determine better the range of 
applicability of velocity centroids, including MVCs. 

We have seen that velocity centroids are only reliable at low 
Mach numbers. At the same time there are other techniques 
available that work for strongly supersonic turbulence, such as 
VCA and VCS. The techniques are complementary and can be 
used simultaneously. For subsonic turbulence, VCA and VCS 
can be used with higher mass species. Note, that the differ- 
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ent techniques pick up different components of velocity ten- 
sor. This can potentially allow a separation of compressional 
and solenoidal parts of the velocity field (LE03). Determining 
these components are crucial for understanding properties of 
interstellar turbulence, its sources and its dissipation. 

We stress that every technique has its own advantages. Ve- 
locity centroids can reproduce velocity better than VCA and 
VCS when the velocity statistics is not a straight power-law. 
Therefore velocity centroids can better pick up dissipation and 
injection of energy scales. 

Nonetheless, one have to bear in mind that spectral indices 
derived from all of the techniques above do not provide a com- 
plete description of turbulence. Anisotropy is another param- 
eter that can be studied. For instance, we showed that velocity 
centroids are useful to study the anisotropy of the turbulent 
cascade, even for highly supersonic turbulence. We showed 
however, that this is restricted to a relatively large mean field. 
For super-Alfvenic turbulence, a model favored by some re- 
searchers (Padoanetal. 2004a), the anisotropy of centroids 
is marginal, which allows to test these theories. Anisotropy 
is not only present in velocity centroid maps, but in other 
statistics as well, su ch as spectral correlation function (see 
iLazarian et aT]l20()l . Combining the two measure can im- 
prove the confidence in the results. 

Additionally, we need tools to study the turbulence vari- 
ations in space (i.e. intermittency). Higher order velocity 
structure functions (those described in this paper are second 
ordert have been shown to be a promising for such studies 
(seelMuller & Biskamp 2000; Cho et al. 2003). Velocity cen- 
troids can be easily recasted in terms of higher order statistics, 
opening a new window for intermittency studies. Is worth 
noting that higher moments can provide the directions of the 
intermittent structures if those get oriented in respect to mag- 
netic field. Our results indicate that such studies can be car- 
ried out even for high Mach number turbulence. In addition, 
studies of intermittency with centroids can be incorporated to 
other techniques. For instance, the interpretation of results 
of the Principal Component Analysis (P CA) tech nique sug- 
gested as a tool for turbulence studies by iHever & Schloerb 
(11997) depends on the degree of intermittency of turbulence 
((Brunt et al. 2003; Heyer 2004 private communication). 

Velocity centroids have been used for studies of turbu- 
lence for many decades. How can we comment on these 
studies from the point of view of our present theoretical un- 
derstanding? Let us glance at the available studies. Tur- 
bulence in H II regions is usually subsonic. Therefore ve- 
locity centr oids could be probably trustworthy there (see 
lO'Dell & C astaneda 1987). Supersonic turbulence in molecu- 
lar clouds (see Miesch & Ballv 1994) is a field where velocity 
centroids might be in error. The same is probably true for 
H I studies (see Miville-Deschenes et al. 2003a)^ where tur- 
bulence is highly supersonic in cold H I. 

We feel that a more careful analysis of particular conditions 
present in the regions under study is necessary, however. One 
one hand, while the turbulence is supersonic for molecular 
clouds, the cores of molecular clouds a re mostly subsonic (see 
iTafalla et al . 2004: Mvers & Lazarian 1998). For such cores 
velocity centroids provide a reliable tool for obtaining veloc- 
ity statistics, if the resolution of the cores is adequate. On the 

^ An interesting feature found in lMiville-Desclienes et ai\ (see 120033) is 
tliat tlie statistics of velocity centroids is almost identical to the statistics of 
density field. This could be due to the dominance of the 71 term, which can 
be tested. 



Other hand, our results are based on the analysis of isothermal 
numerical simulations. In the presence of substantial density 
contrasts caused by the co-existence of different phases the 
velocity centroids may get unreliable even for subsonic turbu- 
lence. Therefore testing of the necessary criteria discussed in 
this paper may be advantageous not only with the molecular 
and H I data, but also with data obtained for H II regions. We 
also note, that for the present analysis we assumed that the 
emission is optically thin and the emissivity is proportional to 
the first power of density. Discussion of more complex, but 
still observationally valuable cases wi ll be done elsewhere. 
The w ork started in this direction (see ILazarian & PogosvanI 
12001 IS encouragmg. 

Rece nt work on centr oids includes papers by 
Ossenkopf & M ac Low! (|2002) where they noticed that 
centroids poorly reproduce velocity statistics for supersonic 
turbulence. Our work above confirms their finding and also 
establishes a regime (Ms < 2.5) when the results by centroids 
can be reliable. 

A quite different and optimistic conclusion about centroids 
was obtained in Miville-Deschenes et al. (2003b). They used 
Brownian noise artificial data and obtained an excellent corre- 
spondence between the centroids and the underlying velocity. 
From the point of view of our analysis the origin of this corre- 
spondence is in the choice of the mean density level. In these 
simulation in order to make the density positively defined the 
authors were adding a substantial mean density to the fluctuat- 
ing density. It is clear from eq.(l29> that adding mean density 
increases of the contribution of /2(R) term that is the term 
responsible for the velocity contribution. 

This work is complementary to the work on delta variance 
of centr oids that we do i n collaboration with Ossenkopf and 
Stutzki dOssenkopf et alJl2005) . There Brownian noise simu- 
lations are used, but extra care is taken to avoid being misled 
by the effect of adding mean density. An iterative procedure 
is proposed there, which allows to correct for the contribu- 
tion of the cross-terms when density-velocity correlations are 
negligible. 

In this paper we are dealing mostly with unnormalized cen- 
troids. The modified centroids (MVCs) allow a different out- 
look at the problem of obtaining velocity statistics from the 
small scale asymptotics. Indeed, our analysis in the paper 
shows that the term /4(R) is unimportant for the cases that we 
studied. In addition, we believe that velocity has steep spec- 
tra. Therefore, if the density is also steep, asymptotically the 
term /3(R), which is then steeper than both velocity and den- 
sity, should be negligible (see Appendix E. 1). As the result, // 
we see from the integrated intensity maps that density indeed 
is steep, we can use MVCs as suggested in LE03, that for suffi- 
ciently small lags are bound to represent the velocity statistics. 
The critical scale at which this is true can be obtained using 
the analytical expressions found in Appendix B.l. Formally 
to find such critical scale one should know the mean density. 
However, if the inertial range is sufficiently long, one should 
not be worried about the exact value for that critical scale. 

We have not seen much advantages of such a use in our nu- 
merical runs because of the limited inertial range available. 
However, unlike numerical simulations, astrophysical condi- 
tions pro vide us with a substantially large r inertial range. For 
instance. lStanimirovic & Lazarianl ( i2001h showed that the tur- 
bulent spectrum for velocity spans from the size of the SMC, 
which is '-^ 4 kpc, to the minimal scale resolved, i.e. ^ 40 pc. 
In partially ionized gas we expect this spectrum to proceed 
to sub-parsec scales. In fully ionized gas (see also discussion 
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in iLazarian et al.lf2004l) the scale of Alfven and slow mode 
damping may be hundreds of kilometers only. 

The limitations of this asymptotic approach stem from the 
fact that for supersonic turbulence the density field tends to 
become shallower. But, as we can see from the MHD sim- 
ulations used here, at A4s ^ 2.5 we find that it is still steep 
and therefore the MVCs should be reliable asymptotically. 
Therefore, for handling observational data we can provide an- 
other prescription: "If the inertial range is sufficiently long 
and density is steep then MVCs asymptotically represent the 
velocity field." For instance, steep density was reported in 
iMiville-Deschenes et al. (2003b). Thus for such density fields 
asymptotic use of MVCs should be advantageous. 

In the cases when the underlying statistics is not a power 
law the interpretation of centroids is more involved. This 
we have seen with our analysis of non-power-law data from 
MHD simulations. While centroids do represent injection 
and damping scales, the integration along the line of sight 
does interfere with a more fine detail study. In this situations 
one should use the full technique of inversion as proposed 
in Lazarian (1995). However, we expect that in most cases 
astrophysical turbulence is self-similar over large expanses of 
scales and therefore the power-law representation is adequate. 

We have made use of the analytical expressions for unnor- 
malized velocity centroids (LE03). Our work as well as a 
study in Levrier (2004), show that normalization of centroids 
improve them rather marginally, while it makes the analyti- 
cal insight essentially impossible. In this paper we tested that 
major conclusions reached for unnormalized centroids are ap- 
plicable to the normalized ones. 

7. SUMMARY 

In this paper we have provided a systematic study of ve- 
locity centroids as a technique of retrieving velocity statistics 
from observations. We used both results of MHD simula- 
tions as well as "reformed" data-sets which have larger in- 
ertial range. We found: 

1. Centroids of velocity can be successfully used to re- 
trieve the scaling properties of the underlying 3D veloc- 
ity field for subsonic turbulence. For supersonic turbu- 
lence with sonic Mach number > 2.5 velocity centroids 
failed to trace the spectral index of velocity. 

2. Our numerics confirmed the expression for centroid 
statistics obtained in LE03. And, in particular, we 
tested successfully the criterion in LE03 for the reliable 
use of velocity centroids. We showed that it reflects a 



necessary condition for centroids to reproduce the ve- 
locity statistics. 

3. We studied modified velocity centroids (MVCs) pro- 
posed in LE03. It is shown shown that MVCs reflect 
the statistics of velocity better than unnormalized cen- 
troids for small lags. This result is valid for both steep 
and shallow density fields with steep velocity. 

4. We showed that velocity-density cross-correlations are 
marginally important for our data set, at least for small 
lags. Combined with the fact that products of density 
and velocity structure functions get subdominant for 
steep velocity and density at small lags, this provides 
a way to reliably study turbulence using MVCs. 

5. We demonstrated that velocity centroids can be used 
for both subsonic and supersonic turbulence to study 
the anisotropy introduced by the magnetic field. Even 
when they fail to retrieve the velocity spectral index in 
supersonic turbulence. For up to at least A^ , ^ 7 they 
for sub-Alfvenic turbulence provide reliably the direc- 
tion of the component of the magnetic field perpendic- 
ular to the LOS. 

6. If turbulence is super- Alfvenic it results in a marginal 
anisotropy of velocity centroids which provides a good 
way to test whether the Alfven Max number of turbu- 
lence in molecular clouds is less or larger than unity. 

7. Within their domain of applicability, centroids of ve- 
locity provide a good tool for turbulence studies that 
should be used in conjunction with other tools, e.g. 
VCA and VCS. 
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APPENDIX 

A. TURBULENCE STATISTICS IN THREE DIMENSIONS 

A. 1 . Turbulence statistics in real (xyz) space. 

The two-point correlation function of a vector field u(x) = [mj:(x), My(x), M2(x)] is defined as: 

Bijir) = {ui(xi)uj(x2)) , (Al) 
where r = X2 -xi is the separation or "lag", in Cartesian coordinates j = x, y, z, and (...) denote ensemble average over all 
space. An additional definition can be made in terms of the fluctuations of the field. This can be obtained formally by replacing 
u(r) by u(x) = u(x)- (u(x)) in equation dAH . Note that this necessarily implies (u) = 0. We will refer to this variation simply as 
correlation function of fluctuations, and denote it by 

Bijir) = (m;(xi)m/x2)) = Bijir) - (m,-) {uj) . (A2) 
In the same manner, it is customary to define the structure function of the same vector field u(x) as 

Dij(r) = ( [m,(xi) - M,(X2)] [m/xi) - m/x2)] ) 

= ([M,(Xl)-Mi(X2)] [m/xi)-m/x2)] ) . (A3) 
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Fij(k) = PND(k) = — ^ / e-'^'-'Bijir) d'^v, (A5) 



Notice that structure functions, by defin ition , dep end only on the fluctuating part of the field and are insensitive to the mean value 
of the field. Combining equations (lAU and ( IA3> it is trivial that: 

Dijir) = 2 [BijiO) - B,./r)] = 2 [B,./0) - Z?,/r)] . (A4) 

Where BijiQ) is also known as the variance. The correlation of fluctuations must vanish at infinity, that is Bijir) ^ 0, as r — > cx). 
Then from equation JA4> is clear that Dijioo) = 2B,j(0). Thus, if we know the structure function of a field we can obtain the 
correlation function of fluctuations and vice versa. However, to get the correlation function in general (as in equation lAll ) we 
also need to know the mean value of the field ((u)). Correlation and structure functions are equivalent in theory. In practice 
however, where we have a restricted averaging space, it is easier to determine Dij{r) more accurately compared to Bij(r). At the 
same time, Bij{0) is usually better determined than D,y[oo] (see discussion in Monin & Yaalom 1975). 

Alternatively, we can use spectral representation to describe turbulence. The spectral density tensor, or N-dimensional power- 
spectrum, Fjj(k) = P^o(k), is given through a Fourier transform of the correlation function of fluctuations; 

1 

(2V)'* 

where k is wave-number-vector. The definitions presented above are general, and also apply to scalar fields. 

The power-spectrum of velocity (in the incompressible regime) has an important physical interpretation as the distribution of 
kinetic energy (per unit mass) as a function of scale. If u(r) is the velocity field, the power-spectrum is the energy in scales 
between k and k+ Sk, and thus the total energy is proportional to (u(r)^) = J PAi/)(k)dk. Note however, that the physical interpre- 
tation is quite different if we talk about the spectra of other quantity (e.g. the power-spectra of density fluctuations). For isotropic 
fields the correlation, structure or spectral tensors can be expressed via longitudinal (parallel to r, denoted by subscript "LL") or 
transverse (normal to r, denoted by subscript 'WA^") components ( Monin & Yaglom. 1975.) : 

(A6) 
(A7) 

(A8) 

where Sij is a Kronecker delta (Sij = 1 for / = j, and Sij = for / j). Solenoidal motions (divergence free, therefore incompress- 
ible) correspond to the transverse components whereas potential motions (curl free, compressible) correspond to the longitudinal 
components. 

A. 2. Turbulence statistics as observed (position-position-velocity space) 

Spectroscopic observations do not provide the distribution of gas in real space coordinates (x = [ji:,3;,z]). Instead we observe 
the intensity of emission of a given spectral line at a position X in the sky, and at a given velocity v along the LOS (we use capital 
letters for 2D vectors and lower case for 3D vectors). Observational data are usually arranged in matrices with coordinates (X, v), 
also called position-position-velocity (or simply PPV) cubes. We will identify the LOS with the coordinate z- Thus, in the plane 
parallel approximation the relation between real space and PPV space is that of a map (X,z) (X, vj, with X = ix,y). 

At any point the LOS velocity can be decomposed in a regular flow, a thermal, and a turbulent components [v^ix) = reg(x) + 
Vtherniai + Vz.turbix)]. This Way, the distribution of the Doppler shifted atoms follows is a Maxwellian of the form: 

where (3 = kbT /m, kb is the Maxwell-Boltzmann constant, T the temperature, and m the atomic mass. In PPV space, the density 
of emitters p.5(X, v,) can be obtained integrating along the LOS 



Bij(r) = 


[Bidr)- 


-BMNir)]"^ 
r'- 


+Bm(r)Sij, 


Dij{r) = 


[DLLir) 


-DMA ^ 


+DNN(r)Sij, 


Fij(k) = 


[FLdky 




+ FNN(k)Sij, 



p,(X,v,)dX dv;: 



dz p(x) (/)(x) 



dX dv„ (A 10) 



where p(x) is the mass density of the gas in spatial coordinates. The d ensity of emitters p^X, v,) can be identified as the column 
density per velocity interval, commonly referred as d^V /dv. Equation jAim simply counts the number of atoms at a position in 
the plane of the sky X, with a z component of velocity in the range [v,, v^ + dv,], and the limits of integration are defined by the 
LOS extent. The integrated intensity of the emission (integrated along the velocity coordinate) corresponds to the column density 
under the assumptions of optically thin media and emissivity linearly proportional to the density: 

/(X) = Ja p,(X, v.) dv, = j a p(x) dz. (Al 1) 

B. PROJECTION OF STRUCTURE FUNCTIONS OF A POWER-LAW SPECTRUM FIELD 

In this section we exemplify the long-wave dominated (steep) case as if coming from a velocity field, while the shallow case 
as if from a density field. The results are interchangeable, depending on the specific spectral index they have (although there is 
no physical motivation to consider a shallow the velocity field). 
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B. 1 . Projection of a field with a steep power spectrum 

Consider a homogeneous, and isotropic velocity field with a steep power-law spectrum. The LOS (chosen to correspond with 
the z direction) velocity structure function will be of the form: 

(K(xi)-y,(x2)]2)=Ci r". (Bl) 



Where r = \/R^ + (z2-Z\)^, R is the separation in the plane of the sky (R^ = [x2-xi]^+ [y2-y\f), and (Z2-Z\), the separation 
along the LOS. This way we can rewrite the 3D power-law structure function as 

{lv-Axi)-v-M2)f) = C,[R' + (z2-z,)T^\ 
The projection of the velocity field, as described in §§2.2 results in 



(B2) 



{mXi)-V2(X2)f) = Ci J~"" {[R^ + (Z2-Zlf]'"'^-(Z2-Zir} dZ2dZi. 



(B3) 



Where Ztot is the largest scale along the LOS. Since the field is isotropic, it is possib le to evaluate this integral changing variables 
from zi and Z2, to z+ = {z\ +Z2)/2 and Z- = Z2-Zi ■ With this new variables, eauation lB3l becomes 



{[V,iXi)-V2(X2)f)=2Q 



-iz-/2) 



= 2Cl 



Z-ll 

(R'+zl] 



{R'+z 



2\m/2 



dz+dz_ 



;i/2 



(z,„,-z-) dz_ 



(B4) 



The remaining integral can be solved analytically in terms of hypergeometric functions, and converge only for m > -1 (which is 
autom atically satisfied since the spectrum is steep therefore with m > 0), yielding (a similar formula can be found in lStutzki et alJ 
[1991: 



1 m 3 
2'~2"'2' 



{[V,(Xi)-V2(X2)f)=2 Ci zi, { R"' 2F1 

zToAm+l)+{zl, + RT^' 



~R^ 



^tot 



R" 



izl+R') 



11/2 



m+2 



(B5) 



(B6) 



m + 2 I 

Where 2F\ is the hypergeometric function, with a series expansion (hypergeometric series) of the form: 

ab X a(a+\)b{b+\)x^ 

2Fx{a,b,c\x) = y{x)=l + — — + — -— ir+ ■■■ 

c 1! c(c+l) 2! 

cy0,-l,-2,-3,... 

However, for all practical purposes, we can calculate it numerically, directly from the definition in equation JB4> . yet the zero point 
(determined by C\) has to be estimated. At small separations {R <^ Ztot), the projected structure functions is well approximated 
by a power-law of the form 

([K(Xi)-y2(X2)]2)«c;/;"'^' R^z,o,. (B7) 

Where Cj is a constant that can be related to C\ by matching the zero point of the 2D structure function from the data to a 
numerical computation using equation (IB3> . For large separations (R 3> Ztot) equation ( IB5> . will follow a power law of the form 



([y,(Xi)-y2(X2)]') « c>"' R » z,o 



(B8) 



WTiere Cj is another constant that will depend on the zero point of velocity fluctuations and on Ziot as well. If we have enough 
inertial range below Ztot- The spectral index m can be obtained from the 2D structure function, and verified with the power 
spectrum. An example of a numerically integrated structure function using equation (IB4> and the asymptotic in equations JB7> 
and ( IB8> is shown in Figure lBTI for a structure function index of m = 2/3 . This first pa nel is almost equivalent to Figure (|3}, where 
we calculate the structure function of integrated data-cubes (Gaussian), but in Figure iBTI onlv theoretical expressions have been 
used, agreement of the two results provides us with a healthy verification. 

B.2. Projection of a field with a shallow power spectrum 
Consider a shallow density field (also isotropic), with a 3D structure function of the form (combining eqs. IIA4I and JSJ) 



((p(Xi)-p(X2))2)=2 



B(0)-C2 



r,/2 



7)12 



The structure function of integrated density (column density) for this case can be written as 



{[I(Xi)-I{X2)Y)=2C2Zto 



r,/2 



v/2 



dz, 



(B9) 



(BIO) 
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Fig. B 1 . — Integrated structure functions from a field with a power-law spectrum. a)The solid line is the exact result (numerically integrating eq.|B3)) for a 
field with steep spectrum ( 3D po wer spectrum index of 730 = —11 /3), the dashed is line for a shallow spectrum (3D power spectrum index 730 = -2.5, obtained 
numerically integrating eq. lBlOl ). The dash-dotted, are the power-law asymptotics for the steep spectrum, and the dotted line for the shallow spectrum. Figure 
b) is the same as a) but we integrated for separations beyond the thickness of the object (chosen here to be 216 to compare with the Gaussian Cubes in Figurel3l. 
we can see clearly the effect of limited thickness discussed in the text. 

where for a shallow spectrum 7/ < 0. This integral is more difficult to evaluate analytically than the long-wave dominated case, 
but for practicality we can solve it numerically. We also find that for small scales, the result can be well approximated by a simple 
power-law: 

([/(Xi)-/(X2)]2)«C;/?'"^' R^Zro,. (Bll) 

Similarly, the zero point can be found by matching the numerical result of eq. (IBl 1> to the calculated (2D) structure function 
from the data. For large separations {R ^ ztot) the resulting structure function of the integrated field will follow the asymptotic 
scaling of the 3D structure function: 



([/(Xi)-/(X2)]^) « constant = c'jR^ R > z,o 



(B12) 



In Figure IBll we also sh ow a n ex ample of the numerical calculation of an integrated structure function using ( IB10> and the 
asymptotics in equations ( IBIH and ( IB12t . for an index rj = -0.5 (corresponding to a 3D power spectrum with an index 73B = -2.5). 

C. POWER SPECTRA OF A FIELD INTEGRATED ALONG THE LINE OF SIGHT 

The procedure presented here can be repeated analogously also for a scalar field (e.g. density). Consider the correlation 
function a velocity, projected along the LOS direction (chosen here to be z): 



(CI) 
(C2) 

(C3) 

(C4) 

(C5) 

(C6) 

(C7) 



(y,(R)y,(X-FR)) = (^y t/zv,(x) j dzv,{^ + {r) ^ 

= jj dzidz2{vz(ii)Vzi^ + r)) = JJ dzidziBzzir). 

(y,(R)V^(X + R)) = Jj dzidz2{vz(x)v,(x + r)) = JJ dzidziB^Ar)- 
Or in terms of the underlying 3D spectrum (inverting eauation lA5> 

(y,(R)y,(X + R))= J J dzxdzi 



d'^ e"' '"K2(k) 

The two-dimensional power spectrum can be obtained taking the Fourier transform of the last expression 



^2d,k(K)=^ / t/^Re 



iK R 



dzidzi 



d\ e'^ 'F.,(k) 



And using k = (K,^.), r = (R,Z2-Zi): 



1 



d'-R e 



iK R 



We can interchange the order of integration: 

^2D,v,(K) =^JJ dzxdzi 
After performing the integral over R one has: 

P2D.vX^)= 1 1 dzidzi I d\' [e'Mz^--''>J(K'-K)F,,(K',£) 



JJdzidzi J d\e'^-^^"'^^''-''^F,,iK,k,) | 
J rf2Re'(K'-K)R+<fe-z,)^^^(K>;) I 
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Now integrate over K : 



P2D.vXK)= dzidzi 



Changing variables to z+ = (zi +Z2)/2, Z- = Z2-Z\, using /v,(-k) = Fzz(k), after of integration in z-. 



(C8) 



P2Dy,(^) = ^^ I dz^ I dk,6{k,)F,,(K,k^) . (C9) 
Now we do the integral over k,, then finally over z+: 

P2D,v, (K) = 4nJ dz^ [F,AK, 0)] 

= 47rz,„,K,(K,0), (CIO) 

lastly, replacing k^ = 0, in equation (IA8> F^z(K,Q) = Fnn{K), we recover from the integrated structure function the spectrum of the 
solenoidal component of the velocity: 

P2D,Vz(^) = '^^Z,o,F^^(K). (Cll) 

D. SECOND ORDER STRUCTURE FUNCTION OF UNNORMALIZED VELOCITY CENTROIDS 
Consider the structure function of the unnormalized velocity centroids (eauation l2H : 

{[SiXi)-S(X2)f) = ^ («/ ^z(xi) P(xi) dz-a J v,(x2) p(x2) dz^ ^ . (Dl) 
As described in §§2.2 (and in lLazarianll995h . we can rewrite last line as in equations J22> and J23> . that we rewrite here: 

{[S(Xi)-S(X2)f)=a^ JJ dzidz2 [Z)(r)-D(r)|x^^xJ - (D2) 

where 

Dir) = ((v;(xi)p(xi)- v,(x2)p(x2))2) . (D3) 
Using the definitions v. = vq + v,, p = po + pin jD3t we have 

Dir) = { [VOPO + Vz(Xl)po + V'op(Xl) + V,(Xi)p(Xi) - VqPO + V,(X2)po + VoP(X2) + V%(X2)p(X2)]^) 

= ({vo[p(xi)-p(x2)] + po[v^(xi)-Vj(x2)] + [v^(xi)p(xi)-v^,(xi)p(x2)]}^^ . (D4) 
Expansion of the last expression yields 

D(r) = {vl [p(Xi) - p(X2)]^ + pI [Vz(Xi) - Vz(X2)]^ + [v,(Xi)p(Xi) - V;:(X2)p(X2)] 

+2voPo [p(Xl) - p(X2)] [Vz(Xi) - V,(X2)] + 2po [Vz(Xi) - Vz(X2)] [v,(Xi)p(Xi) - V^(X2)piX2)] 
+2vo [v,(Xi) - V,(X2)] [Vj(Xi)p(Xi) - Vj(X2)p(X2)]) 
= Vo ([p(Xl)-p(X2)]^) +/30 ([v,(Xi)- V,(X2)]^) + ([02(Xi)p(Xi)-{^,(X2)p(X2)]^) + 

2vopO ([p(Xl) - p(X2)] [Vz(Xi) - V2(X2)]) + 2po {[v^iXi) - V^(X2)] [v,(Xi)p(Xi) - Vj.(X2)p(X2)]) 

+2vo ([p(Xi) - p(X2)] [v,(Xi)p(Xi) - V,(X2)p(X2)]) . (D5) 

Now consider the cross-terms one by one, the third term in equation ( ID5I I is 

{[v,(xi)i5(xi)-v,(x2mx2)f) = {v^Axi)p^^^^^^ (D6) 

Relating the fourth order moments with the second order using the Milhonshikov hypothesis ^ (see iMonin & Yaelomll 19751) 

{hih2h-ihii) « {h\h2) {hT,h/i) + {hih-i) {hjh/i) + {hxhji) {h2h^) we have 

(v2(xi)p2(xi)) = (v,(xi)p(xi)y,(xi)p(xi)) « {tixi)) {p\xi)) + 2{v,{xi)p{x,)f , (D7) 

(v2(x2)p2(x2)) = (v,(X2)p(X2)v,(X2)p(X2)) « (v^^^^)^ (p2(X2)) +2 (y,(X2)p(X2))' , (D8) 
2 (Vz(Xi)p(Xi)v,(X2)p(X2)) W 2(v,(Xi)p(Xi)) (Vz(X2)p2(X2)) 

+2(v,(Xi)Vj(X2))(p(Xi)p(X2)) 

+2(y,(xi)/5(x2)) (v,(x2)p(xi)) . (D9) 

* Evidently, this liypotliesis is identically true for Gaussian fields, strongly non Gaussian fields may show deviations. 
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If we combine the last three lines using 



( [Vz(Xi) - V,(X2)f) = (vliXi)) - 2 (v,(Xi)v,(X2)) + (V'2(X2)) 

= 2(v%2)-2(y,(xi)v,(x2)), 
(y,(xi)y,(x2)) = (v2) - i ([v,(xi)-y,(x2)]2) , 

(p(Xi)p(X2)) = (p')- ^ ([p(Xi)-p(X2)]2) . 



(DIO) 
(DID 

(D12) 



Together with {v,(xi)p(xt)) = {vz(x2)p(x2)) =0, {v^(xi)) = {v^ixi)) = (v?), and (p^(xi)) = (p^ixz)) = {f?) equation (ID6> reduces 



to 



([0,(Xi)p(Xi)-0,(X2)p(X2)]2)«2(v2)(p2)-2(v,(Xi)/5(X2))(v,(X2)p(Xi)) 



(v/)--([v,(xi)-y,(x2)]2) 



(p2)-_([p(xi)-p(X2)]2) 



(D13) 



To treat terms of the form {pv,) we can generalize the correlation function defined in equation JAIL considering a four- 
dimensional vector of the form [v^(x), v^.(x), Vz(x), p(x)]. The resulting cross-correlation between the z component of the velocity 
and the density is 

B,p(r)=(v,(xi)p(x2)). (D14) 

Similarly to the derivation of equation \A6\ it can be shown that for an isotropic (four-dimensional) field Bjp(r) = BLp(r)rj/r, and 
Bpj(r) = Bpi{r)rj/r = -Bip{r)rj/ r. Thus {v,{x\)p{x2)) = -{v.(x2)p(xi)), and equation (ID13> simplifies to 

([v,(Xi)p(Xi)-y,(X2)p(X2)]')«2(y,(Xi)p(X2))2 + (vT2)([p(xi)-p(X2)]') 

+ (p')(K(xi)-y,(x2)]') 

([v.(Xi)-V,(X2)]') {[p(Xi)-p(X2)f) . 



The next term of equation {D5i 

2voPo ([p(Xl) - p(X2)] [Vj(Xi) - V,(X2)]) = 0, 

as shown explicitly in lMonin & Yaglomljl975l) . The fifth term in equation SD5\ is: 

2po ([Vz(Xl)- Vz(X2)] [Vj.(Xi)p(Xi)-y.(X2)p(X2)]) =2po {v^(Xi)p(Xi)) + 2po (v.(X2)p(X2)) 

-2po {v^(xi)v^(x2)p(xi)) 

-2po {v^(Xi)v^(X2)p(X2)) 

For the terms with correlations of third order we need to introduce the so called two-point third order moments 

Bijjir) = (m,(xi)m;(xi)m,(x2)). 



(D15) 
(D16) 



(D17) 



(D18) 



Which can be generalized to include cross-correlations of density and velocity as explained for equation (ID14> . Analogously to 
the second order cross-correlations, considering an isotropic (four-dimensional) field, and decomposing it in terms of longitudinal 
and lateral components. It can be shown (for more details we refer the reader to Monin & Yaglom 1975) that 



(i'z(Xl)Vz(X2)p(Xi)) = (v,(Xi)v,(X2)p(X2)) , 
(>'z(Xl)Pz(Xl)p(X2)) =- (Vz(Xi)P2(Xi)p(X2)) . 

Equation ( ID19> , combined with {p(xi)v^{xi)) = (/5(x2)v?(x2)) reduces equation (ID17> to: 

2po ([vz(xi)- v,(x2)] [v,(xi)p(xi)- Vj.(x2)p(x2)]) =4/90 (p(xi)y^(xi)) 

-4/90 (v,(xi)v,(x2)p(xi)) . 

Similarly the last term in equation ( ID5> can be written 

2v() {[p(xi)-pix2)] [v,ixi)p(xi)-v,(x2)p(x2)]) =2vo (v^(xi)p^ixi)) 

+2vo(vz(X2)/9^(X2)) 
-2vo (/0(Xi)p(X2)i>z(Xi)) 
-2vo (/5(Xi)p(X2)Vz(X2)) . 

But in this case , (v.(xi)p^(xi)) = (vz(x2)p^(x2)) = 0, combined with equation ( ID20> , yields 

2vo ([p(Xi)-p(X2)] [Vz(Xi)p(Xi)-i5,(X2)/5(X2)]) = 0. 



(D19) 
(D20) 



(D21) 



(D22) 
(D23) 
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Finally, combining (ID5IID15l . lDT6l . lD2niD23t 

D(r) « vl { [p(xi) - p(x2)f) + pI { [v,(xi) - v,(x2)]') + 2 (v,(xi)p(x2))' 

+ {vl) {[p(Xi)-p(X2)f) + {p^) ([v,(Xi)-V,(X2)]') 

-\ (K(Xl)-Vz(X2)]') {[p(Xi)-p(X2)f)+4po{p(Xt)vl(Xi)) 

-4po (v,(xi)v,(x2)p(xi)) . (D24) 

Grouping some terms: 

D(r)^{p^ + pl) {[v,(xi)-v,(x2)f) + {v^ + vl) {[p(xi)-p(x2)f) 
-\ ([v.(xi)-y,(x2)]2) ([p(xi)-p(x2)]') + 2(v,(xi)p(x2))' 

+Apa {p{xi)vl{xi)) -4po (v,(xi)v,(x2)p(xi)) . (D25) 
And lastly, with (p^) = {fp' + p^), and (v^) = (v' + Vq) we arrive to equations i24\ and J26> : 

D(r)^{p^){[v,(xi)-v,(x2)f) + {v^){[p(xi)-p(x2)f) 

-\ ([v.(xi)-v,(x2)]') ([p(xi)-p(x2)]')+c(r), (D26) 

with 

c'(r) = 2 (v,(xi)p(x2))' + 4p„ (p(xi)v^(xi)) -4po (v,(xi)v,(x2)p(xi)) . (D27) 

We should note that the s econ d term in equation ( ID27> has no contribution to the structure function of centroids because when 
substituted into equation JD2> cancels ((p(xi)v,(xi)^) = (p(xi)v,(xi)^) |xi=X2)- Therefore we omitted this term in the body of the 
paper. Moreover, this term does not appear in the correlation function of unnormalized centroids, not including it in the definition 
of c(r) shows better the symmetry between structure functions with correlation functions and spectra. 

E. CENTROIDS AND THE "CROSS-TERM" FOR POWER-LAW STATISTICS 

We have already discussed the differences between power-law long-wave and short-wave dominated fields {steep and shallow 
spectra). In particular, the fact that although in all cases power spectra are power-laws, but structure functions are only power- 
laws for steep spectra whereas correlation functions are only power-laws (at small scales) for shallow spectra. In this section we 
will sketch the expectations for the cross-term (73 [R]) for the idealized case of infinite power-law spectra, considering all the 
different shallow and steep combinations. 

E. 1 . Steep density and steep velocity 

In this case the structure function of both fields is well described by a power law. Consider power-law isotropic underlying 
statistics of the form (i,.,(r) = C^r'", dp{r) = Cpr'\ where m, n > 0, and C,,, Cp are constants. Here the cross-term will be 

/3(R) = -^a^QCp J I AziAz2 [r'"^" - r'""" lx,=xj : (El) 

which is analogous to the projection of a steep field with a spectral index 73/5 = -(m + n)-3. The resulting cross-term is negative 
and steeper than both density and velocity. Thus at sufficiently small scales its contribution can be safely neglected compared to 
the integrated density and velocity. 

E.2. Shallow density and steep velocity 

Here the structure function of velocity is a power-law d^Xr) cx C^r™, with m > and C,, constant. For small separations 
(relative to a critical scale as discussed in the main body of the paper) the density fluctuations will have a power-law correlation 
function (p(xi)p(x2)) ~ Cpr", where n < and Cp is constant. Therefore the structure function of density can be presented as 
dp{r) K, 2({fp-) -Cpr"). For this combination of indices the cross-term becomes 

/3(R) «-(p2) {[V(Xi)-V(X2)f) + a^Q.CpJ j Az,Az2 [r™^"- r'«+"lx,=xj • (E2) 

The first term on the right, although negative has the same slope as /2(R). The two terms inside the integrals are equivalent to 
the projection of a field with an spectral index 730 = -(m + n)-3. Which in this case is going to be shallower than the velocity, 
and could be either positive or negative. Posi tive for m+n > 0, and negative for m+n < 0. If the level of density fluctuations (p^) 
is large the first term in the right of eq. (IE2t will dominate the slope of /3(R). However such term will be canceled with part of 
/2(R), changing effectively the weighting of the integrated velocity from (p^) = Pq + (p^) to Pq. The resulting structure function 
of unnormalized centroids can be written as^ 

' Equations IE3IIE5I and lE7l coincide witli the decomposition we derived in lOssenkopf et all 120051) . 
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( [5(Xi) - S(X2)f) ^{vl+{vl)){ [/(Xi) - /(X2)]') + pI { [V (Xi) - y (X2)]') 

+a^Q.Cp f[dzidz2 r'«+"lx:=xj +'f4(R). (E3) 



The exact contribution of the terms inside the double integrals to /3(R), and in general the importance of the cross-term to the 
centroid statistics will depend on the constants C,, and Cp. If at large scales the magnitude of the term containing the structure 
function of integrated velocity and the cross-terms are comparable, then at small scales the centroids could fail tracing velocity. 

E.3. Steep Density and shallow velocity 

This case is somewhat analogous to that described above. Consider for the density a power-law structure function dp(r) cx 
Cpr", with n > and Cp constant. For the velocity, a power-law correlation function that translates into a structure function 
t/,, (r) w 2((v^) -Ci.r'") with m <Q and C,. constant. We have already discussed that a shallow velocity is not physically motivated. 
However if we have a steep velocity field without infinite power-law behavior we could have a structure function that resembles 
the shallow case. For instance if the structure function "saturates" at large separations it will be similar to that of a shallow field, 
which grows rapidly at small scales and then flattens at large scales. The corresponding cross-term is 

/3(R) « -(0?) {[I(Xi)-I(X2)f) + a^Q.CpJJ dzidzi [r'"^"- r"'^"lx,=xj • (E4) 

It has a component that scales as the column density (first term on the right). The remaining component could be positive (if 
m + n > 0) or negative (if m+n < 0). Unnormalized centroids for these indices can be expressed as 

([5(Xi)-5(X2)]2)«vg([/(Xi)-/(X2)]2) + (p^+(/52))([y(Xi)-y(X2)]2) 

WQ.Cp J I dzidz2 [r'""" - r'"- I^^^^J + ^(R). (E5) 

E.4. Shallow density and shallow velocity 

For this combination of indices we can consider power-law correlation functions (for small separations). Yielding structure 
functions of the form dp{r) « 2((/5^) - Cpv"), dvAr) ~ 2((v'^) -Cyr"^), with m, n < Q, and Cp, Cy constants. As a result the 
cross-term becomes 

/3(R) ~-{vl){ [/(Xi) - /(X2)]') -{p^){ [V(Xi) - V(X2)f) 

-la^CCpj I dzidzi ['•'"^"- '•'"^"lx,=xj ■ (E6) 

Now we have a term that scales as the column density, a term that scales the integrated velocity, and the term inside the integrals 
(now negative). Notice also that for large values of the velocity or density dispersion, the contribution of column density or 
integrated velocity respectively is increased within /3(R). But such contributions will be canceled for the unnormahzed centroids 
that will have a structure function: 

( [5(Xi) - S{X2)f-) « vl { [/(Xi) - I(X2)f) + pI { [V(Xi) - V(X2)f) 

-la^CyCp J I dzidzi [r'"^" - r'"^" j^^^xj +^4(R). (E7) 

F. CORRELATION FUNCTION AND POWER SPECTRA OF UNNORMALIZED CENTROIDS 

Structure functions are considered to be preferable statistics according to lMonin & Yaglomljl975h . They are subjected to less 
errors related to averaging. In any case, correlation functions are trivially related to each other by the formula in equ ation( IA4> . 
Using the expression for the structure function of unnormalized centroids from LE03 one gets (see also lLevrieil2004l for a direct 
derivation): 

{S{Xi)S{X2))^{p'){vl)(az,o,f + a^vl jj {p(xi)p(x2))dzidz2 + a^pl jj (v,(xi)v,(x2)) dzidzj 
jj {pixi)p(x2)) (V2(xi)v2(x2)) dzidz2-a^ jj {p(iii)vzix2)fdzidz2 
+2 a^po j j [(p(xi)v,(xi)0(x2))] dzidz2. (Fl) 

In analogy with equation (I27> we can rewrite ( IF1> as: 

(5(Xi)5(X2)) « (p2)(v2)(az,„,)2 + Bl(R) + B2(R) + B3(R)+M(R), (F2) 
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where 



BliR) = a^vl jj (p(xi)p(x2))dzidz2, 

B2(R) = aVo jj (vz(xi)v,(x2))dzidz2, 

B3iR) = a^ jj (p(xi)p(x2)) (v\(xi)y,(x2))dzidz2, 

B4{R) = -a^jj {p{xi)v,(x2)fdzidz2 + 2 a^o jj [(p(xi)v,(xi)v(x2))] dzidzz- 



(F3) 
(F4) 
(F5) 
(F6) 



Here B1(R) is Vy times the correlation function of fluctuations of column density, B2(R) is a^p^ the correlation function 
integrated velocity, B3(R) is different from 2/3(R) only by a constant, and B4(R) contains the density-velocity cross-correlations 
found in /4(R). 

The Power spectrum of unnormalized centroids is the Fourier Transform of equation iFl\ : 

P2D.s(^) = {p^){vl)(az,orfS{K) + vlP2Dj(K) + a^plP2DyX^) + T{B3(R)} + T{B4(R)}. (F7) 

It is interesting of this is that the weighting factors are no longer (v^ + v^) and {p\ + p-), but only v\ and p\. And in fact the 
contribution of the column density alone can be eliminated by removing the mean LOS velocity vq. However the cross-term, and 
density-velocity cross-correlations do affect the scaling properties of the centroids of velocity. And as the turbulence increases in 
strength and the ratio {p') / p\ grows, the "contamination" will also be stronger (see Ossenkopf et al. 2005). 

The correlation function of MVCs can be obtained by subtracting (v?) times the correlation function of column density from 
the correlation function of normalized centroids. Or equivalently: 

CFmvc{^) = (5(Xi)5(X2))- (/(Xi)/(X2)) 

« {p^){vl)(az,o,f + B2(R) + B3\R) + B4(R), (F8) 

with 



B3'(R) = -a^jj (pixiM^D) [{vl)-{v,(xiWx2))]dzidz2 
■ jj {pixi)pix2))dyXr)dzidz2. 



a 



(F9) 



And, thus the power spectrum of MVCs is: 

P2D.Mvc(^) = (p'> {vl)(azr„rf6(K) + a Vo^2D,v,(K) + t[b3'{R)} + T {B4{R)} . (FIO) 

As the three measures (spectra, correlation, and structure functions) are trivially related, they represent the same statistics. 
They can be used interchangeably as dictated by practical convenience. 

G. UNNORMALIZED CENTROIDS VS. MVCS 

When are MVCs better than regular unnormalized centroids? To answer this question let us consider MVCs in terms of 
correlation functions. First, we can minimize the contribution of density fluctuations to the unnormalized centroids by setting 
vo = 0. In this case we can fold the difference of unnormalized and modified centroids in the cross-terms B3{R), and B3 (R) 
respectively. Thus, the criterion for MVCs to be an improvement over unnormalized centroids would be |B3(R)| > \B3 R)|. In 
other words: 

I// (p(Xl)p(X2)) (V;(Xi)v,(X2))dzidZ2| 



// {p{xi)p(x2))dyXr)dzidz2\ 



> 1. 



(Gl) 



Let us c onsider a more realistic steep velocity field, with a power-law structure function that saturates at some cutoff scale r,, o 
(see LPOO. lLazarian & PogosvaiJ2004 : 



(G2) 



We use r,, o to avoid confusion with the critical scale for the short-wave dominated spectra, and because this cutoff scale can also 
be related to a correlation length for the correlation function: 



(v,(Xi)v,(X2)) = (v2) 



v,0 



(G3) 



If we combine this steep velocity field, with a steep density field -i.e. (/5(xip(xi)) = (p^)'"p (>/('"" + '"p())-, with the same change 
of variables as before (z_ = Z2-Z1) the criterion in equation ( IG1> becomes 



rZtoi 

Jo 


\ 1 






dz_ 






(jj2+j2)™/2+r;;.^ 


CZtot 

Jo 






(/f2+z2)"/2 


dz- 






(«2+z2)'"/2+r;"„ 



> 1. 



(G4) 
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c) m = 0.1, n = 0.1 
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Fig. G1. — Criterion for MVCs to trace velocity better tlian unnormalized centroids for an example of steep density and velocity, (see equation IG4H The 
parameters are labeled in the title and axis for each panel, 7.tot = 100 in all cases. The criterion is met for parameters that correspond to iso-contours above 1. 
Panels (a), {h), and (c) show the low sensitivity of the position of the 1 iso-contour with the spectral indices. Panels (d), (e), and (f) show the strong dependence 
on r,, o and the lag Ry, as well as a weak dependence on rp q. 



In Figure IgTI we plotted iso-contours of this ratio for various combinations of parameters. Within the parameter space corre- 
sponding to iso-contour above 1.0 MVCs are expected to trace the statisti cs of velocity better than unnormalized centroids. In 
the first three panels (a, b, c) we present, for rp o = r,,,o how the ratio in ea.( IG4> varies with the lag R over a extreme contrast of 
possible spectral indices (within the long-wave dominated regime). In the last three panels (d, e, f) we choose a particular set 
of spectral indices and relax the condition rp o = A,,o. In all cases z,oi = 100. We find that the criterion for MVCs to be better 
than unnormalized centroids is satisfied for lags smaller than the cutoff for the velocity field, with very little sensitivity to the 
particular spectral index or the cutoff scale of density. 

We performed the same exercise for the other physically motivated case, that is, steep velocity and shallow density. A shallow 
density can be approximated by the same expression for the correlation function of density, but with a small rp o. We tested for 
spectral indices from 0.1 to 1.0 for the velocity, and from -1.0 to -0.1 for the density. The same behavior as for the steep-steep 
case was found. That is, MVC are expected to trace the statistics of velocity better than unnormaUzed centroids for lags smaller 
than the velocity correlation-scale (r,. o), and smaller than the LOS extent of the object (ztot)- 
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